Log In Sign Up

Multiple changepoint detection for periodic autoregressive models with an application to river flow analysis

In river flow analysis and forecasting there are some key elements to consider in order to obtain reliable results. For example, seasonality is often accounted for in statistical models because climatic oscillations occurring every year have an obvious impact on river flow. Further sources of alteration could be caused by changes in reservoir management, instrumentation or even unexpected shifts in climatic conditions. When these changes are ignored the statistical results can be strongly misleading. This paper develops an automatic procedure to estimate number and locations of changepoints in Periodic AutoRegressive models. These latter have been extensively used for modelling seasonality in hydrology, climatology, economics and electrical engineering, but there are very few papers devoted also to changepoints detection, moreover being limited to changes in mean or variance. In our proposal we allow the model structure as a whole to change, and estimation is performed by optimizing an objective function derived from the Information Criterion using Genetic Algorithms. The proposed methodology is brought out through the example of three river flows, for which we built models with possible changepoints and evaluated their forecasting accuracy by means of Root Mean Square Error, Mean Absolute Error and Mean Absolute Percentage Error. The last years of data sets have been omitted from the selection and estimation procedure and were then used to forecast. Comparisons with literature on river flow forecasting confirms the efficiency of our proposal.


page 1

page 2

page 3

page 4


Regularization methods for the short-term forecasting of the Italian electric load

The problem of forecasting the whole 24 profile of the Italian electric ...

Stream-Flow Forecasting of Small Rivers Based on LSTM

Stream-flow forecasting for small rivers has always been of great import...

Operator-theoretic framework for forecasting nonlinear time series with kernel analog techniques

Kernel analog forecasting (KAF), alternatively known as kernel principal...

TILDE-Q: A Transformation Invariant Loss Function for Time-Series Forecasting

Time-series forecasting has caught increasing attention in the AI resear...

AdaEnsemble Learning Approach for Metro Passenger Flow Forecasting

Accurate and timely metro passenger flow forecasting is critical for the...

1 Introduction

Discontinuities are often introduced into climatic or hydrologic time series as a result of anthropogenic impacts or changes in instrumentation or location. Further plausible reasons are modifications in reservoir system management or new water pricing. As defined by LL12 a changepoint is ”a time where the structural pattern of a time series first shifts”. In many cases, changepoints are located at known times (dam construction, measuring instrument change) and it is easy to take into account their effects. When changepoints are located at unknown times and their features are ignored, the time series estimation can be misleading (LL07; LWLRGF07). In fact, an undetected changepoint can lead to: misinterpretation of the model, biased estimates and less accurate forecasting (Hansen01). Sometimes, even an abrupt shift in climatic conditions does not have a significant impact on hydrologic time series (BCRCSAS15). Taking all these into account, changepoint detection becomes a demanding job especially if its identification is required soon after occurrence (e.g. flood predictions).

For the past four decades several techniques have been employed for changepoint detection. As far as hydrological applications are concerned, many authors have considered the problem of detecting a single changepoint (Cobb78; Buishand84; HM94; RT96), but very few analyzed more realistic multiple changepoints situations. For example, SLX10 proposed a segmented regression with constraint method in order to model both trend analysis and abrupt change detection. ZLMWRT15 studied changing trends and regime shifts of streamflow using long term historical records at different hydrological stations in the Yellow River basin. A Bayesian method is used to investigate the possibility of changes in the mean values and to locate the changepoint in a hydrological time series  (RT96). Kehagias04 and KNP06

adopt hidden Markov model, and dynamic programming algorithm, respectively, to locate multiple changepoints in hydrological and environmental time series. A nonparametric estimation of both number and positions of changepoints in multivariate settings is proposed by 

MJ14. KS12

use a nonparametric algorithm to solve the problem of changepoint detection by comparing probability distributions over two consecutive intervals.

SLW17 and GSL17 use, respectively, a dynamic factor model and dynamic programming for studying the segmentation of multivariate time series. In order to assess the dependence of hydrological drought variables, SS10 construct a multivariate probability distribution using copulas and estimate its parameters using genetic algorithms. The cumulative sum (CUSUM) approach has proven to be an efficient method to identify the change of annual maximum rainfalls in southern Taiwan (CPL12). Several methods (Pettitt, CUSUM, Segment Neighbourhood) for identifying multiple changepoints are presented by PG17; they also show how the presence of autocorrelation may affect these methods. The problem of modeling a class of nonstationary time series using piecewise Autoregressive processes is considered in DLRY06. A more general procedure for piecewise threshold autoregressive models (TAR) is developed by YTL15. For a recent review of changepoint analysis in time series, see AH13.

Models accounting for seasonality, as the seasonal autoregressive integrated moving average (SARIMA), developed originally by BJ70, have been broadly used in the literature of hydrologic models (MD05; Du10; WGCZ14). Periodic time series models have been introduced because these latter cannot be filtered to achieve second-order stationarity, and this is because the correlation structure of these time series depends on the season (Ve85a; Ve85b; Mc93). MBF17 also showed that the seasonal differencing maintains seasonal correlation structure, whereas periodic term is completely removed by seasonal standardization or by spectral analysis. Using seasonal differenced models can deteriorate the forecasting ability of river flows (DTK76) which is a main problem for planning and operating of water resources. Also extensions of periodic models have been employed in hydrological applications (Per17). General overviews of periodic models and their applications is presented in HM94FP04MW06.

This paper proposes a Periodic AutoRegressive (PAR) model for time series recorded monthly with multiple changepoints at unknown times, which is well-tailored to hydrologic and climatic analysis as the majority of time series observed in these areas are periodic with respect to time. Changepoints detection procedures for such series have been studied by LWLRGF07; LLL10 among others. However, PAR models analyzed in these research papers allow mean changes, while the model studied here allows also PAR and trend parameters to switch at each changepoint time. Our results share a number of similarities with their finding allowing at the same time a generalization of results. Moreover, the portmanteau test (Mc94) designed for diagnosing the accuracy of PAR models is used for each segment. We employ a procedure based on Genetic Algorithms (GAs) to detect changepoints and estimate resulting PAR models, allowing also to identify more parsimonious subset selections. These kind of methods are well suited for complex global optimization, as they have been widely applied to hardly tractable identification and estimation problems (SS10; UP15). General proposals of changepoint detection for time series by means of GAs can be found in JK13; DFHW17 among others.

Our method is employed to analyze the average monthly flows of three rivers: the Garonne (France), and two Canadian rivers, South Saskatchewan (measured at Saskatoon) and Saugeen (measured at Walkerton). Forecasting ability of resulting models is then evaluated by standard performance measures. We found that models with changepoints (possibly due to both climatic reasons and human intervention) outperforms standard PAR models in terms of prediction accuracy. This could led to improvements in operating the reservoir system, which can result in saving large amount of money (an overview of optimization methods used in reservoir operation is presented in FEJ13). The rest of the article is organized as follows: Section 2 describes proposed methodology for estimation and forecasting; data and discussion on the results are included in Section 3; comments close the paper in Section 4.

2 Methodology

2.1 Model description and estimation

We consider the problem of modeling a non-stationary time series by segmenting the series into several PAR processes linked at different changepoints. The period of series is and is assumed to be known. Observation in season of year is denoted by , with and .

There are different segments, each of which contains an integer number of years, and denotes the first year of segment . The first segment includes years from to , second segment contains years from to , third segment contains years from to , and so on. If we indicate with the number of changepoint times, the segment structure is defined as follows:

In order to ensure reasonable estimates, it is required that each segment contains at least a minimum number of years, therefore for any segment . We let , so that if year belong to set then time is in segment . For the sake of simplicity we assume that the total number of observations is a multiple of .

The model driving our work is given by:


where and

The process is a PAR given by:


We assume that trend parameters and depend only on the segment, whereas means are allowed to change also with seasons. The autoregressive maximum model order at season in the -th segment is given by , so that , , represent the PAR coefficients during season of the -th segment. For simplicity, we assume that maximum autoregressive order . In our procedure coefficients will be allowed to be set to zero, in order to get more parsimonious models. This is accomplished by introducing in each segment

a binary vector

of length , named PAR lags indicator, which specify presence or absence of parameters. The first digits represent the subset PAR model for period , subsequent digits are related to period and so on.

The error process in equation (2

) is a periodic white noise, with

and , , , Unless otherwise stated we assume that each segment is periodic stationary with period , in the sense that

for all integers and . Periodic stationarity is discussed in Gl61.

The number of changepoints , the changepoints location and PAR lags indicator are named as structural parameters. They can take discrete values and their combinations amount to a very large number. Once such parameters are determined, the trend intercepts , the slopes , the seasonal means , the AR parameters and the innovation variances (for segment , season and lag ) are analytically estimated. These parameters, assuming that structural parameters are known, can be estimated according to the following steps:

  1. Trend parameter vectors and

    are estimates by Ordinary Least Squares (OLS) method:

    that leads to detrended data

  2. Seasonal means are computed on resulting detrended data as follows:

    and implies: .

  3. AR parameters are estimated separately for each segment and season. Each of these specific series is selected from and is incorporated in a design matrix of dimensions , which includes lagged observations. Parameter constraints are specified by a matrix , where is the number of free parameters. These constraints are designated on the basis of PAR lags indicator as follows:

    • For each lag , the element of vector is evaluated

    • If value is equal to 1 then a row equal to the -th row of identity matrix is added to .

    Final estimate of is obtained by constrained optimization, with linear constraint given by . Explicitly (in matrix form):

    where is obtained by OLS estimation.

  4. Lastly, estimation of innovation variances is performed for each segment and season on final residuals, considering that each regime has a possibly different sample size.

Selection of optimal structural parameters, on the other side, is a complex problem for which no closed form solution is available. As far as it involves the evaluation of a very large number of possible combination, GAs are naturally suitable for this issue.

2.2 Identification of structural parameters

The GA is a nature-inspired optimization method, often employed when it is required to find an optimal solution from a prohibitively large discrete set. In a maximization problem it allows to approximate the optimal solution, which maximizes an objective function (named fitness in GA terminology), by a simple procedure. In the generic iteration a population of binary encoded solutions (called chromosomes) is subdued to genetic operators: the selection randomly chooses chromosomes for the subsequent steps, usually proportionally to their fitness value; by crossover two solutions are allowed to combine together, with a fixed rate , exchanging part of their values and creating two new individuals; lastly, the mutation step allows each binary value to flip its value from 0 to 1, or vice versa, with a fixed probability , providing a further exploration of search space. The resulting population replaces the previous, and the flow of generations stops if a certain condition is met, for example a fixed number of generations. It is also possible, adopting the elitist strategy, to maintain the best chromosomes found up to the current iteration, irrespective of operators effect.

Two different genetic approaches have been analyzed in pilot experiments for our structural parameters identification problem. The first method is a hybrid GA (HGA), in which number of changepoints and changepoint locations are encoded in a binary chromosome, while PAR lags indicator are obtained by enumerating all possible subset models, given each segment and season, and returning only the best one. This latter task is computationally feasible when maximum autoregressive order is small, as far as models must be evaluated for each segment and season. In this way of proceeding optimal values of structural parameters are obtained combining an exact method (exhaustive enumeration) with an approximation procedure (GA). Second algorithm is a simple GA (SGA) in which subset model indicators are included in chromosome. In this case the size of search space and the length of chromosome considerably increase with . As a result of this it was necessary to increase the number of iterations for SGA. Both of these algorithms, however, seemed to be equivalent on exploring the solution space, and for small values of they very often reach an equivalent solution. As far as in this paper we shall consider values of up to , because autoregressive procedure must capture the short term dependence, we employed HGA in our analysis. For larger values of we recommend the use of SGA.

In HGA binary chromosomes encode a candidate segmentation as follows: first two bits give number of changepoints (limited to a maximum of 3 in our study, so that a number of segments up to 4 is allowed); subsequent bit intervals, whose length is custom fixed, produce changepoint times . This part of encoding must ensure following constraints:

due to the fact that a minimum number of observations must be contained in each segment. In order to accomplish this the bit intervals encode real numbers , constructed to determine percentage of remaining values to be attributed to the corresponding -th segment. In fact, when placing a new changepoint there are some illegal positions, due to above specified constraints: this implies that observations must be left out from both the beginning and the end of considered segment. This strategy depends on candidate number of segments, so changepoints are uniquely identified in four possible ways:

  • If (one segment) then .

  • If (two segments) then

  • If (three segments) then:

  • If (four segments) then:

Such an encoding procedure, introduced in BP12, allows each possible chromosome to be legal, so there is no computational time wasted on evaluating infeasible solutions.

As far as fitness function measures goodness of solutions in GAs, our model identification problem shall include a term linked to the goodness of fit and a part related to a penalization on number of parameters. Many options are available: we shall consider a criterion inspired by NAIC, introduced by To90 for threshold models, given by:


where is the model residual variance of series in segment and season , is related to the sample size, is related to the number of parameters, is the penalization term. We adopt a scaled exponential transformation of as fitness function, for maximization purposes: , where is a problem dependent constant which is introduced for suitably scaling fitness.

Before leaving the estimation procedure, we can make a comment. A two-step procedure for identification and estimation of PAR models based on Yule-Walker approach was presented by MW06

. We propose an automatic procedure for detecting changepoints, identifying the best model and estimating the parameters for each segment. Our estimation method is based on multiple linear regression. Both estimation techniques are efficient from statistical point of view.

2.3 Forecasting and performance assessment

The forecasting method employed is the standard one-step-ahead procedure. We shall remove last year from dataset (corresponding to 12 observations) in order to estimate model parameters, and use those data points for evaluating forecasting performance. The logarithm of data is adopted as Box-Cox transformation, performed before fitting the model, as far as it is the most used in monthly river flow analysis and it ensures that model residuals are approximately normal distributed and homoscedastic 

(EV92; MG13).

Forecasting accuracy of proposed models is evaluated with respect to the following measures:

where and are, respectively, true and predicted values of last year observations, on logarithmic scale. These measures are explicitly defined in HK05. A smaller value indicates a better model performance. The MAE and the RMSE between the proposed model and the observed data are calculated in the same units of the observed data. The MAPE measure is based on percentage errors; these latter have the advantage of being scale-independent, and they are frequently used to compare forecast performance between different data sets. These criteria must be interpreted only as an indication that shows which model performs better, but no statement can be made from this comparison. All measures are calculated using the hydroGOF and/or forecast packages in R software. For other measures of forecast accuracy we refer to HK05; KBB05.

2.4 Model validation

The analysis of the residuals for the estimated model is needed to test its relevance. The stationarity of the residuals allows us to apply the standard 95% confidence limits (that is ). To test the joint statistical significance of the residual autocorrelations, Mc94 proposed a Ljung-Box portmanteau test for PAR models:


where designs the period, is the autocorrelation coefficient with lag , is considered the maximum time lag and where represents the integer part of the real number

. The test statistic

follows approximatively a chi-square distribution

, with degrees of freedom.

3 Results and discussion

We shall now study effectiveness of proposed methodology in river flow analysis. Data related to three rivers having different lengths, means of annual flows and located in different regions, will be examined. They consist of:

  • flows of Garonne river measured at Tonneins, France;

  • flows of South Saskatchewan river measured at Saskatoon, Canada;

  • flows of Saugeen river measured at Walkerton, Canada.

Mean, standard deviation and correlation functions of these rivers vary significantly from month to month. Further details will be given in subsequent subsections.

Garonne river has been analyzed in UP15. The South Saskatchewan and Saugeen river flows series are discussed in NMH85, and they are available at They used nine different models to generate thirty-six one-step-ahead forecasts for the logarithmic flows. In terms of RMSE, PAR models give the best results. We will show that our methodology can improve their results.

Estimation and forecasting will be performed as discussed in subsection 2.2 and 2.3. As far as choice of genetic operators is concerned we propose standard roulette wheel selection, bit-flip mutation, and a modified single-point crossover: instead of allowing all bits in chromosome to be selected as possible cutting points, we shall consider only bits that subdivide segmentation . In such a way parameter structures can be naturally inherited by offspring, avoiding the possibility of destroying solutions. Elitist strategy is also employed.

In fitness function, related to NAIC criterion (3), we fixed so that it resembles the penalization structure of AIC. Several experiments have been conducted considering various combinations of parameters , and , in order to provide a variety of explicative models. We used the following strategy: conditioning on four possible values of , which include the model with no changepoints and situations with possible structural breaks up to, respectively, , and , we selected four models for which the best value of fitness function has been observed. Forecasting accuracy of these models, labelled as PAR, will be then evaluated.

Computations will be performed by using Matlab and R softwares.

3.1 Garonne river

The Garonne river which flows through Spain and France is the third largest river in France in terms of flow. Its total length is about 647 km with a catchment area of 51500 km at Tonneins. It is the main contributor to the Gironde Estuary which is the major European fluvial-estuarine system. Flow measures are recorded at the Tonneins gauging station, where there is no tidal effect. Data are obtained from daily discharge measurements in cubic meter per second (m/s) from January 1959 to December 2010 (DIREN-Banque Hydro, French water monitoring). Daily data flows are then transformed in monthly data consisting in flows averaged for one month. The data and the log transformed data (624 observation or 52 years) are presented in Figure 1.

Figure 1: Monthly flows (up) and logarithmic monthly flows (down) for the Garonne river.

Optimal PAR models and forecasting results are reported in Table 1. The first proposed model is a PAR without changepoints which is extensively used to perform hydrological forecasting (NMH85; MW06). All the other proposed models allow one or more changepoints.

Figure 2: Changepoint detected on years 1977 and 1989 for Garonne river

Year 1989 is detected as possible changepoint in almost all configurations. According to CVMHNLLB07, years 1988-1989 seem to be the driest in the 1980-1990 decade. Moreover, the air temperature over Western Europe showed an abrupt shift at the end of the 1980s. For a better understanding of climatic changes and their impact on water resources, BCRCSAS15 studied a subset of 119 temperatures, 122 rainfalls and 30 hydrometric stations over the entire France. They detected a shift in annual mean air temperature in 1987-1988 for more than 75% of the 119 temperature stations. They also detected a shift between 1985 and 1990 for 18 hydrometric stations.

A related series (Dordogne river flows) is available in order to help draw conclusions. Using flow series for Dordogne river from 1959 to 2010, we found exactly the same 2 changepoints (1977 and 1989) as with PAR model for Garonne river. The Dordogne and Garonne rivers are the main rivers of the Adour-Garonne river basin. The climate in the basin area is under the influence of oceanic conditions in its western parts and we suppose that both rivers experience similar climatic conditions. Changepoint corresponding to 1977 could be linked to local effects. Generally, it is difficult to attribute any of the changepoint time to a specific climatic event or anthropogenic factor.

Years of changepoint
PAR / 0.2476 0.2215 2.3562 1.2022
PAR 1989 0.2731 0.2133 2.2665 1.2113
PAR 1977, 1989 0.2725 0.2131 2.2643 1.2146
PAR 1970, 1988, 1998 0.3149 0.2511 2.6100 1.2241
Table 1: Results of evaluation criteria of the logarithmic forecast errors for Garonne

In terms of goodness of fit we observe that fitness values are monotonically increasing with number of changepoints; forecasting accuracy, on the other side, is maximum in terms of MAE and MAPE for model with 2 changepoints corresponding to 1977 and 1989, while best value of RMSE is observed for model with no changepoints.

Figure 2 graphically displays the segmentation with 2 changepoints and Figure 3 shows the logarithmic flows and their forecasts using a PAR model with two changepoints. From latter figure we can see that the RMSE value could be influenced by large error prediction of first observation.

Figure 3: Logarithmic flows of Garonne (full line) and one-step best PAR forecasts (dashed line).

As a diagnostic check, the residual autocorrelations up to lag 36 were computed. Figure 4 shows the autocorrelation function (ACF) for the residuals of PAR model with 2 changepoints. The two full lines indicate lower and upper bounds of the ACF assuming that residuals are white noise. The residuals for PAR model with 2 changepoints indicate no correlation. This graphic checking provides ”some evidence” on the adequacy of proposed PAR model. Validity of our model is also evaluated with the portmanteau test given in equation (4) and results are presented in Table 2. All P-values in this table suggest that the proposed model is not rejected at the usual 5% significance level, except for October in the first segment, with an empirical significance level of about 2.5%. However, this does not strongly point to the model inadequacy as that empirical level is still not too low and, consequently, the possibility of an error of type I exists.

Figure 4: Autocorrelation function (ACF) of the residuals of the fitted PAR model with two changepoints to the Garonne flow.
1959-1977 1978-189 1989-1999
January 0.4217 0.7622 0.6915
February 0.1701 0.2600 0.5702
March 0.3938 0.5050 0.1058
April 0.9323 0.2758 0.7183
May 0.6557 0.0789 0.6045
June 0.0761 0.1109 0.5058
July 0.6058 0.7771 0.5916
August 0.3658 0.6296 0.1500
September 0.3581 0.2201 0.8310
October 0.0251 0.0950 0.7068
November 0.3318 0.1811 0.6034
December 0.1558 0.3261 0.5973
Table 2: P-values of the portmanteau test defined in eq. (4) with .

3.2 South Saskatchewan river

The South Saskatchewan River originates in the Rocky Mountains and drains an area of about 139,600 km. It joins the North Saskatchewan river and form Saskatchewan river, which is the fourth longest river in North America.

The time series of mean monthly flows of the Saskatchewan river measured at Saskatoon, Canada from January 1912 to December 1976 include 780 observation (65 years). The monthly series and logarithmic monthly flows can be viewed in Figure 5.

Figure 5: Monthly flows (up) and logarithmic monthly flows (down) for the South Saskatchewan river.

NMH85 fitted the transformed data with nine models and they recommended the PAR model. Using a PAR without changepoints we obtain a similar model as in NMH85; the slight differences arise from the fact that our model allows intermediate constraints and we use more data. In addition, our method was applied to estimate a PAR model with at least one changepoint. Optimal PAR models are showed in Table 3.

Figure 6: Changepoint detected on years 1937, 1961 and 1968 for Saskatchewan river

In all models with at least one possible changepoint we detect year 1968. January 1969 corresponds to a modification in the reservoir system management: the Gardiner Dam came into full operation upstream from Saskatoon on the South Saskatchewan River. The Gardiner Dam is the third largest embankment dam in Canada and one of the largest embankment dams in the world. The reservoir provided valuable benefits to the community: power generation, recreational benefits, and also magnitude of floods have been lessened, with minimum flows downstream guaranteed throughout the year. Before the creation of Gardiner Dam, snowmelt in conjunction with summer rains led to heavy flooding.

To model the effect of the operation of the Gardiner dam on the average monthly flows of the South Saskatchewan river HMM77 developed an intervention model. To see how the monthly average changes after the intervention, the cumulative sum was calculated and plotted for each month. Using flows measured at Saskatoon from 1942 to 1974, they found that the operation of the Gardiner dam significantly affected the average monthly flows: the average flows increase from November to March and decrease from April to September (Figure 7).

Figure 7: Monthly means of the South Saskatchewan river before (full line) and after (dashed line) of the construction of Gardiner dam

Looking into the history, the changepoint corresponding to 1937 could be linked to the drought conditions on the Saskatchewan prairies during the Dust Bowl years of the 1930s. The drought came in three waves, 1934, 1936, and 1939-1940, but some regions experienced drought conditions for as many as eight years. The changepoint corresponding to 1961 could arise as an effect of the beginning of dam construction in 1964.

As for the Garonne river, the performance of models was considered in terms of fitness, as far as goodness of fit is concerned, and RMSE, MAE and MAPE for evaluating accuracy of one step-ahead logarithmic forecasts (Table 3). The PAR model with 4 segments (PAR) performs better in terms of both estimation and forecasting in comparison with other models. Figure 6 display the logarithmically transformed series with 3 detected changepoints, and Figure 8 shows the logarithmic flows and their forecasts using such model.

Years of changepoint
PAR / 0.66428 0.5126 10.2718 1.2253
PAR 1968 0.4360 0.3167 6.3646 1.2553
PAR 1937, 1968 0.4348 0.3115 6.2757 1.2643
PAR 1937, 1961, 1968 0.4275 0.2933 5.8488 1.2661
Table 3: Results of evaluation criteria of the logarithmic forecast errors for Saskatchewan

To check for the whiteness of the residuals the autocorrelations up to lag 36 were computed for the PAR model. Only one significant autocorrelation was found (Figure 9). Table 4 shows P-values of portmanteau test: they suggest that the proposed model is not rejected at 5% significance level, except for August 1961-1968, November 1938-1961 and December 1969-1975 where the P-values are about 1.5%, which does not strongly suggest model inadequacy. In this situation, PAR model without changepoints seems inappropriate at the 5% nominal level.

Figure 8: Logarithmic flows of South Saskatcewan (full line) and one-step PAR forecasts (dashed line).
Figure 9: Autocorrelation function (ACF) of the residuals of the fitted PAR models with three changepoints to the South Saskatchewan flow.
1912-1937 1938-1961 1961-1968 1969-1975
January 0.9536 0.6893 0.5412 0.0423
February 0.3981 0.5351 0.3869 0.2332
March 0.8564 0.5745 0.1022 0.5783
April 0.3221 0.0622 0.4431 0.4106
May 0.6451 0.5700 0.0985 0.1154
June 0.5260 0.9274 0.6738 0.7060
July 0.5384 0.5004 0.3646 0.5631
August 0.1553 0.3217 0.0180 0.0885
September 0.1081 0.3099 0.4733 0.3288
October 0.4983 0.8576 0.5301 0.4271
November 0.1029 0.0137 0.1615 0.4746
December 0.7882 0.0841 0.4347 0.0144
Table 4: P-values of the portmanteau test defined in eq. (4) with .

3.3 Saugeen river

The Saugeen River is located in southern Ontario, Canada; it begins in the Osprey Wetland Conservation Lands and flows generally north-west about 160 kilometres (99 miles) before exiting into Lake Huron. Starting from 1950 it is served by Saugeen Valley Conservation Authority (SVCA), a corporate body founded for managing and preserving water and other natural resources in river watershed. It is also among the most important areas in Ontario for fishing, specifically nearby significant dykes as Denny’s Dam, in which a popular conservation area is located.

Data analyzed are average monthly flows from January 1915 until December 1976, measured at Walkerton, Ontario.

Figure 10: Monthly flows (up) and logarithmic monthly flows for Saugeen River
Model Years of changepoint
PAR / 0.4851 0.3712 11.1849 1.1871
PAR 1970 0.3520 0.2869 9.0176 1.1918
PAR 1965, 1970 0.3756 0.2961 9.3380 1.1944
PAR 1950, 1958, 1970 0.3762 0.2966 9.2641 1.1972
Table 5: Results of evaluation criteria of the logarithmic forecast errors for Saugeen

Table 5 shows results of optimal PAR models with changepoints: year 1970 is always detected as changepoint in our models. One possible reason would be related to works aiming to reconstructing Denny’s Dam. In fact, between the end of 1960s and the beginning of 1970s, Great Lakes Fishery Commission managed to rebuild Denny’s Dam in order to provide an effective bloackage against parasites such as sea lamprey, preventing them from infiltrating in river. Being Denny’s Dam among the biggest dykes of river course this could have been a non-ignorable effect on its flow. There have also been important human work on Saugeen in the 1950s: the main reason of SVCA creation in 1950 was, indeed, flood control management. Walkerton business district, which is the gauging station, has been subject to major floods in early and mid 1900. This has led to the construction, starting from 1956, of 2.4 km of dykes and floodwalls to protect the central business district as well as residential neighborhoods from potential floods.

As far as goodness of fit is concerned, we observe that also in this case fitness is monotonically increasing with number of changepoints; in terms of forecasting the best performance is observed for model with single changepoint in 1970 considering all measures. Figure 11 plots time series with this changepoint.

As far as this time series has already been modeled in literature we can also perform some comparisons. WIZX07 proposed a functional-coefficient autoregression model in order to estimate and forecast monthly flows of Saugeen. Forecasting performance, measured on natural data, have been compared to PAR(1) model by HM94, resulting in an improvement in terms of from to . Corresponding value of our model computed on natural data is , which further improves standard PAR(1) model.

Autocorrelation function (ACF) for the residuals from the PAR model with one changepoints is reported in Figure 12 and it shows very few significant autocorrelations. P-values of portmanteau test for the model are reported in Table 2. These values indicate that model is not rejected at significance level, except for September of first segment.

1915-1969 1970-1976
January 0.6072 0.7949
February 0.7640 0.2823
March 0.9022 0.4626
April 0.5054 0.3544
May 0.1796 0.2453
June 0.8815 0.6919
July 0.4122 0.6855
August 0.5044 0.2895
September 0.0129 0.2006
October 0.3257 0.3474
November 0.8046 0.3445
December 0.6125 0.1522
Table 6: P-values of the portmanteau test defined in eq. (4) with .
Figure 11: Changepoint detected on year 1970 for Saugeen
Figure 12: Autocorrelation function (ACF) of the residuals of the fitted PAR models with three changepoints to the Saugeen flow.

4 Conclusions

The goal of our research was to develop a computational procedure to estimate the number of changepoints and their locations in time series with periodic structure. We analyzed river flow series, for which accurate forecasting is a highly demanding issues in hydrology, especially for the management of reservoir systems. An undetected changepoint can lead to much less accurate forecasting. Hence identifying a changepoint becomes an important task as the rapid environmental change act differently on river flows: the magnitude and timing of river flows change, and new patterns of extreme droughts or floods are observed.

Our procedure estimated changepoints for time series related to river flows of Garonne (France), Saugeen and South Saskatchewan (Canada). The reasons for such changepoints are possibly due to both human activities and climatic oscillations. The comparison of our findings with forecasting results of NMH85 and WIZX07 confirms efficiency of the proposed method.

In our study, we examined monthly data with changepoints allowed only at the end of the year (that is, a multiple of number of seasons). Modifications of method proposed in the present paper are under study: techniques for monthly, weekly or daily time series with periodic structure allowing changepoints at any season are worth pursuing. Therefore, detecting a changepoint in the middle of a year will prevent dispersing its effects over adjacent seasons. Moreover, as far as PAR models are based on a large number of parameters, one could question on whether it is necessary to consider a separate AR model for each season: we allowed to build subset PAR models in order to conveniently decrease number of parameters, but a considerable gain in parsimony would be achieved by reducing number of seasons in PAR model (FP04 and HM94

proposes several statistical hypothesis tests).

Lastly, it is known that a stationary autoregressive processes has a short memory (BD13; Ro03). Time series which exhibit long range dependence are characterised by autocorrelations which decay very slowly, while a stationary autoregressive process have rapidly decaying autocorrelations. Hydrological data or internet traffic data generally exhibit structural changes and long range dependence (SB13). Therefore a long memory process with periodic structure could be appropriate for hydrological data.

The authors thank professor Francesco Battaglia for his valuable and constructive remarks. Part of this work has been supported by the project Green Water Labex Cote.

Conflict of Interest Statement
The authors declare that they have no conflict of interest.