# An Application of Filtered Renewal Processes in Hydrology

Authors: Mario Lefebvre; Fatima Bensalma
Journal: International Journal of Engineering Mathematics
 (2014)

---

## Abstract

Filtered renewal processes are used to forecast daily river flows. For these processes, contrary to filtered Poisson processes, the time between consecutive events is not necessarily exponentially distributed, which is more realistic. The model is applied to obtain one- and two-day-ahead forecasts of the flows of the Delaware and Hudson Rivers, both located in the United States. Better results are obtained than with filtered Poisson processes, which are often used to model river flows.


---

## Body

## 1. Introduction

Let{
N
(
t
)
,
t
≥
0
} be a Poisson process with rate λ. A filtered Poisson (sometimes calledshot noise) process is a continuous-time stochastic process {
X
(
t
)
,
t
≥
0
} defined by
(1)
X
(
t
)
=
∑
n
=
1
N
(
t
)
‍
w
(
t
,
τ
n
,
Y
n
)
(
X
(
t
)
=
0
if
N
(
t
)
=
0
)
,
in which the random variables τ
1
,
τ
2
,
… denote the arrival times of the events of the Poisson process and Y
1
,
Y
2
,
… are assumed to be independent and identically distributed random variables that are also independent of {
N
(
t
)
,
t
≥
0
}. The function w
(
·
,
·
,
·
) is called the response function.In many applications, the response function is chosen of the form(2)
w
(
t
,
τ
n
,
Y
n
)
=
Y
n
e
-
(
t
-
τ
n
)
/
c
,
where c is a parameter that must be estimated. It then gives the value at time t of an event of magnitude Y
n of the Poisson process that occurred at time τ
n. Moreover, the random variables Y
1
,
Y
2
,
… are generally assumed to be exponentially distributed with parameter μ. With the above response function, the filtered Poisson process behaves as in Figure 1. Actually, this behavior depends on the form of the response function, but not on the distribution of the time between the events. Therefore, the same behavior would be observed in the case when {
N
(
t
)
,
t
≥
0
} is a renewal process. Remember that a Poisson process is a particular renewal process.Figure 1
Example of a trajectory of a filtered Poisson process defined by (1) and (2).To model the daily flows of rivers given their most recent observed values, conceptual and physical models based on the filtered Poisson process have been widely used successfully for many decades; see, for example, Weiss [1], Kelman [2], Koch [3], and Konecny [4]. Filtered Poisson processes are still being used to model various phenomena in civil engineering; see Yin et al. [5] and Miyamoto et al. [6].Now, especially in hydrological applications, the form of the response function in (2) is taken for granted rather than being justified by the observations of the variable of interest. Actually, it is generally simply an assumption made to obtain a mathematically tractable model.Similarly, the actual distribution of theY
n’s is not investigated. However, if one is only interested in forecasting the value of X
(
t
+
δ
), based on the values of the process up to time t, this does not really cause a problem because the forecast is normally based on the mean of the Y
n’s, and the distribution itself is not needed.Finally, and even more importantly, the assumption that{
N
(
t
)
,
t
≥
0
} is indeed a Poisson process is not tested. Again, this is a simplifying assumption, because one can then make use of the nice properties of the Poisson process, notably the fact that it has independent and stationary increments.Next, notice that, with the response function in (2), the effect on the value of the process of an event that occurred at time τ
n is maximum at that time instant and {
X
(
t
)
,
t
≥
0
} has discontinuities at the arrival times of the events of the Poisson process. In the case of the application in hydrology that consists in modeling the flow of rivers, if one looks at real hydrographs, one sees that there is a more or less extended period of time during which the flow increases from a minimum to a maximum. Even if river flows data are generally available on a daily basis, one does not observe very sudden increases of the flow, followed by an exponential decrease. Rather, it takes on average two to three days for the flow to reach a maximum, and the decrease is more or less rapid.To obtain a more realistic model for the variations of river flows, by taking into account the periods of flow increase, Lefebvre and Guilbault [7] used the following response function:
(3)
w
(
t
,
τ
n
,
Y
n
)
=
Y
n
(
t
-
τ
n
)
d
e
-
(
t
-
τ
n
)
/
c
,
where d is a positive parameter. With this response function, the river flow begins to increase immediately after time τ
n, but the maximum is reached after d
c time units (i.e., generally days) and then it starts to decrease. Their work was improved by the authors, who provided a method to estimate the parameter d. They found that, for the applications that they considered, d was in the interval (0,1).With a Poisson process being a continuous-time Markov chain, the time that the process spends in a given state is exponentially distributed, which is a strictly decreasing density function. LetT
1
,
T
2
,
… be the times between two successive events, that is, the interarrival times. In the case of river flows, the interarrival time is generally a random variable with a density function that is increasing at first. Hence, the assumption that {
N
(
t
)
,
t
≥
0
} is a Poisson process is at least doubtful.In Lefebvre [8], the author tested the hypothesis that the T
n’s are exponentially distributed. He found that, for the Delaware River, it was not accurate. Instead, it was shown that both a gamma and a Rayleigh distribution fitted the data much better. The objective was to forecast the peak flows of the Delaware River, which is a difficult task due to the very weak correlation between the successive peaks. To try to forecast these peak flows, the author used a filtered renewal process as a model for the river flow. An advantage in using a Rayleigh distribution for the interarrival time is that the filtered renewal process can be transformed into a filtered Poisson process by working on a different time scale.Poisson and renewal processes have been used in hydrology to model low flows (see, for instance, Loaiciga and Leipnik [9] and Yagouti et al. [10]) as well as various physical phenomena such as traffic noise (see Marcus [11]). Andersen [12, 13] gave new theoretical results on filtered renewal processes (see also the references therein).The aim of the present paper is to show that we can obtain more accurate short-term forecasts of river flows by making use of filtered renewal processes rather than filtered Poisson processes. To do so, we will first need to develop a formula to forecast the flow at timeX
(
t
+
δ
), given the history of the process up to time t. Because we cannot appeal to the memoryless property of the exponential distribution, the task of finding an estimator for X
(
t
+
δ
) is more difficult. The values of the river flow before time t must be taken into account, contrary to the case of filtered Poisson process, for which the forecasts depend only on X
(
t
).Moreover, because we want to clearly see the improvement obtained by using a more realistic distribution for the interarrival time, we will take the response function defined in (2). Indeed, as mentioned above, filtered Poisson processes having this response function are commonly used in hydrology, and the forecasting formulas are then easy to derive and implement. Therefore, it will be easier to judge the quality of the model that we propose.In the next section, filtered renewal processes will be formally defined and the formulas needed to forecast the river flows will be developed. Then, in Section3, the results will be applied to forecast the flows of the Delaware and Hudson Rivers. Finally, a few concluding remarks will be made in Section 4.


## 2. River Flows Modeled as a Filtered Renewal Process

LetX
(
t
) denote the flow of a certain river at time t. We assume that
(4)
X
(
t
)
=
∑
n
=
1
N
(
t
)
‍
Y
n
e
-
(
t
-
τ
n
)
/
c
(
X
(
t
)
=
0
if
N
(
t
)
=
0
)
,
in which {
N
(
t
)
,
t
≥
0
} is a renewal process. The τ
n’s are the arrival times of the events of the process {
N
(
t
)
,
t
≥
0
}, and the times T
n
=
τ
n
-
τ
n
-
1 between the successive renewals are independent and identically distributed random variables. Finally, the (positive) random variables Y
1
,
Y
2
,
… are independent and identically distributed (and independent of {
N
(
t
)
,
t
≥
0
}).Thus, the process is characterized by the sudden increases of the river flow caused by the events that occur at timesτ
1
,
τ
2
,
…. A renewal cycle is defined as the time T
n between two consecutive peaks. For the interarrival time, the following distributions (in addition to the exponential distribution) will be considered in Section 3: lognormal, gamma, Weibull, and Rayleigh.Next, in order to be able to forecast the river flow at timet
+
δ, we must first estimate the parameter c that appears in the response function. If there are no signals between t and t
+
1, then X
(
t
+
1
) is given by
(5)
X
(
t
+
1
)
=
e
-
1
/
c
X
(
t
)
.Hence, to estimate the parameterc, we can consider all the days where the flow decreased in the data set and calculate the ratio X
(
t
+
1
)
/
X
(
t
). The point estimate of c is then derived from the arithmetic mean R
¯ of this ratio:
(6)
c
^
=
-
1
ln
⁡
R
¯
.In the case of the random variablesY
n, which constitute a random sample of a parent random variable Y, their exact distribution is not needed. Indeed, only the expected value of Y will be used to forecast the value of X
(
t
+
δ
). Therefore, we simply have to compute the mean increase of the river flow during the calibration period of the model.We will now derive formulas to estimate the river flow one and two days in advance. These formulas can be used with any distribution for the interarrival timeT, which denotes the parent variable of the T
n’s.


### 2.1. Forecasting the Flow at Timet
+
1

We will first try to forecast the river flow at timet
+
1, given all the information available up to t. To do so, we assume that there will be at most one event in the interval (
t
,
t
+
1
]. At any rate, as mentioned above, in practice flow values are generally available on a daily basis. Therefore, it is not really possible to determine whether more than one event occurred in (
t
,
t
+
1
]. When it does happen, these events can be considered as a single large event.Apart from the case when the interarrival timeT has an exponential distribution, we must take the history of the process into account to calculate the probability that there will indeed be an event in (
t
,
t
+
1
]. Let
(7)
U
i
=
τ
N
(
t
)
+
i
-
τ
N
(
t
)
+
i
-
1
for
i
=
1,2
,
…
.
That is, U
i is the time elapsed between the (
i
-
1
)
st and the i
th event after the most recent one that occurred before or at time t. The random variables T and U
1 are identically distributed. We must compute
(8)
p
k
:
=
P
[
U
1
≤
k
+
1
∣
U
1
>
k
]
,
where k is the (integer) number of days since the most recent signal was observed. We have
(9)
p
k
=
∫
k
k
+
1
‍
f
U
1
(
u
)
d
u
∫
k
∞
‍
f
U
1
(
u
)
d
u
.
With the help of a mathematical software, we can obtain numerical values of this probability for any density function f
U
1 and any value of k.Next, ifT is exponentially distributed, it is well known that, given that there is a signal (i.e., an event) between t and t
+
1, the conditional density of T is uniform on the interval (
t
,
t
+
1
], so that the signal in question occurred on average at time t
+
1
/
2. In the general case, we must compute the conditional density function of U
1 in the interval (
k
,
k
+
1
], given that U
1
>
k.If there is indeed a signal betweent and t
+
1, the forecast of the river flow at time t
+
1 would be
(10)
Z
k
:
=
X
(
t
)
e
-
1
/
c
^
+
E
[
Y
]
1
P
[
U
1
>
k
]
×
∫
k
k
+
1
‍
e
-
(
k
+
1
-
u
)
/
c
^
f
U
1
(
u
)
d
u
.
Therefore, the forecast of the river flow at time t
+
1, given the entire history of the process up to t, is given by
(11)
X
^
(
t
+
1
)
=
X
(
t
)
e
-
1
/
c
^
(
1
-
p
k
)
+
Z
k
p
k
=
X
(
t
)
e
-
1
/
c
^
+
p
k
E
[
Y
]
1
P
[
U
1
>
k
]
×
e
-
(
k
+
1
)
/
c
^
∫
k
k
+
1
‍
e
u
/
c
^
f
U
1
(
u
)
d
u
.


### 2.2. Forecasting the Flow at Timet
+
2

LetU
1
,
k
=
U
1
∣
{
U
1
>
k
}. To forecast the river flow at time t
+
2, given all the information available up to time t, we need the distribution of the random variable S
k
:
=
U
1
,
k
+
U
2, where U
2 is defined in (7).We now assume that the probability that there will be more than two signals in the interval(
t
,
t
+
2
] is negligible (irrespective of the value of k). The density function of S
k is given by the convolution of the density functions of U
1
,
k and U
2:
(12)
f
S
k
(
s
+
k
)
=
∫
k
∞
‍
f
U
1
,
k
(
u
)
f
U
2
(
s
+
k
-
u
)
d
u
.
Since U
2 is a nonnegative random variable, we can write that
(13)
f
S
k
(
s
+
k
)
=
∫
k
s
+
k
‍
f
U
1
,
k
(
u
)
f
U
2
(
s
+
k
-
u
)
d
u
for
s
>
0
.In the case whenT has an exponential distribution with parameter λ, we know that (
S
k
-
k
) has a gamma distribution with parameters α
=
2 and λ. In the other cases, we must generally evaluate the above integral numerically.The value ofX
(
t
+
2
) is given by(14)
X
(
t
+
2
)
=
{
e
-
2
/
c
X
(
t
)
if
U
1
,
k
>
k
+
2
,
e
-
2
/
c
X
(
t
)
+
Y
1
*
e
-
(
k
+
2
-
U
1
,
k
)
/
c
if
U
1
,
k
≤
k
+
2
and
U
2
>
2
-
(
U
1
,
k
-
k
)
,
e
-
2
/
c
X
(
t
)
+
Y
1
*
e
-
(
k
+
2
-
U
1
,
k
)
/
c
if
U
1
,
k
+
U
2
≤
k
+
2
+
Y
2
*
e
-
(
2
-
U
2
)
/
c
(
and
U
3
>
2
-
(
U
1
,
k
-
k
)
-
U
2
)
,where Y
i
* denotes the size of the i
th signal in the interval (
t
,
t
+
2
]. As mentioned above, we assume that
(15)
P
[
U
3
>
2
-
(
U
1
,
k
-
k
)
-
U
2
]
≃
1
.We must first compute(16)
p
1
,
k
:
=
P
[
U
1
,
k
∈
(
k
,
k
+
2
]
]
.
When U
1 has an exponential distribution with parameter λ,
(17)
p
1
,
k
=
P
[
U
1
∈
(
0,2
]
]
=
1
-
e
-
2
λ
.
In general, we have
(18)
p
1
,
k
=
∫
k
k
+
2
‍
f
U
1
(
u
)
d
u
∫
k
∞
‍
f
U
1
(
u
)
d
u
.Next,(19)
p
2
,
k
:
=
P
[
U
1
,
k
+
U
2
∈
(
k
,
k
+
2
]
]
≡
P
[
S
k
∈
(
k
,
k
+
2
]
]
=
∫
0
2
‍
f
S
k
(
s
+
k
)
d
s
=
∫
0
2
‍
[
∫
k
s
+
k
‍
f
U
1
,
k
(
u
)
f
U
2
(
s
+
k
-
u
)
d
u
]
d
s
.Alternatively, using the fact thatU
1
,
k and U
2
(
≥
0
) are independent random variables, we may write that
(20)
p
2
,
k
:
=
P
[
U
1
,
k
+
U
2
∈
(
k
,
k
+
2
]
]
=
∫
k
k
+
2
‍
P
[
U
1
,
k
+
U
2
∈
(
k
,
k
+
2
]
∣
U
1
,
k
=
u
]
f
U
1
,
k
(
u
)
d
u
=
∫
k
k
+
2
‍
P
[
U
2
∈
(
k
-
u
,
k
+
2
-
u
]
]
f
U
1
,
k
(
u
)
d
u
=
∫
k
k
+
2
‍
[
∫
0
k
+
2
-
u
‍
f
U
2
(
v
)
d
v
]
f
U
1
,
k
(
u
)
d
u
.Finally, by independence,(21)
P
[
{
U
1
,
k
≤
k
+
2
}
∩
{
U
2
>
2
-
(
U
1
,
k
-
k
)
}
]
=
∫
k
k
+
2
‍
P
[
U
2
∈
(
k
+
2
-
u
,
∞
)
]
f
U
1
,
k
(
u
)
d
u
=
∫
k
k
+
2
‍
[
∫
k
+
2
-
u
∞
‍
f
U
2
(
v
)
d
v
]
f
U
1
,
k
(
u
)
d
u
=
∫
k
k
+
2
‍
[
1
-
∫
0
k
+
2
-
u
‍
f
U
2
(
v
)
d
v
]
f
U
1
,
k
(
u
)
d
u
=
p
1
,
k
-
p
2
,
k
.The forecast of the river flow at timet
+
2, given the entire history of the process up to time t, is thus given by
(22)
X
^
(
t
+
2
)
=
X
(
t
)
e
-
2
/
c
^
(
1
-
p
1
,
k
)
+
{
X
(
t
)
e
-
2
/
c
^
+
E
[
Y
]
∫
k
k
+
2
‍
e
-
(
k
+
2
-
u
)
/
c
^
f
U
1
,
k
(
u
)
d
u
}
×
(
p
1
,
k
-
p
2
,
k
)
+
{
X
(
t
)
e
-
2
/
c
^
+
E
[
Y
]
∫
k
k
+
2
‍
e
-
(
k
+
2
-
u
)
/
c
^
f
U
1
,
k
(
u
)
d
u
+
E
[
Y
]
∫
k
k
+
2
‍
[
∫
0
k
+
2
-
u
‍
e
-
(
2
-
v
)
/
c
^
f
U
2
(
v
)
d
v
]
kkkkkkkkkkkkk
×
f
U
1
,
k
(
u
)
d
u
[
∫
0
k
+
2
-
u
‍
e
-
(
2
-
v
)
/
c
^
f
U
2
(
v
)
d
v
]
}
p
2
,
k
=
X
(
t
)
e
-
2
/
c
^
+
{
E
[
Y
]
∫
k
k
+
2
‍
e
-
(
k
+
2
-
u
)
/
c
^
f
U
1
,
k
(
u
)
d
u
}
p
1
,
k
+
{
E
[
Y
]
∫
k
k
+
2
‍
[
∫
0
k
+
2
-
u
‍
e
-
(
2
-
v
)
/
c
^
f
U
2
(
v
)
d
v
]
l
l
l
l
l
l
l
l
l
l
×
f
U
1
,
k
(
u
)
d
u
∫
0
k
+
2
-
u
‍
f
U
2
(
v
)
d
v
}
p
2
,
k
=
X
(
t
)
e
-
2
/
c
^
+
{
E
[
Y
]
e
-
(
k
+
2
)
/
c
^
∫
k
k
+
2
‍
e
u
/
c
^
f
U
1
,
k
(
u
)
d
u
}
p
1
,
k
+
{
E
[
Y
]
∫
k
k
+
2
‍
[
e
-
2
/
c
^
∫
0
k
+
2
-
u
‍
e
v
/
c
^
f
U
2
(
v
)
d
v
]
k
k
k
k
k
k
k
k
k
k
k
k
×
f
U
1
,
k
(
u
)
d
u
∫
0
k
+
2
-
u
‍
e
v
/
c
^
f
U
2
(
v
)
d
v
}
p
2
,
k
.


## 2.1. Forecasting the Flow at Timet
+
1

We will first try to forecast the river flow at timet
+
1, given all the information available up to t. To do so, we assume that there will be at most one event in the interval (
t
,
t
+
1
]. At any rate, as mentioned above, in practice flow values are generally available on a daily basis. Therefore, it is not really possible to determine whether more than one event occurred in (
t
,
t
+
1
]. When it does happen, these events can be considered as a single large event.Apart from the case when the interarrival timeT has an exponential distribution, we must take the history of the process into account to calculate the probability that there will indeed be an event in (
t
,
t
+
1
]. Let
(7)
U
i
=
τ
N
(
t
)
+
i
-
τ
N
(
t
)
+
i
-
1
for
i
=
1,2
,
…
.
That is, U
i is the time elapsed between the (
i
-
1
)
st and the i
th event after the most recent one that occurred before or at time t. The random variables T and U
1 are identically distributed. We must compute
(8)
p
k
:
=
P
[
U
1
≤
k
+
1
∣
U
1
>
k
]
,
where k is the (integer) number of days since the most recent signal was observed. We have
(9)
p
k
=
∫
k
k
+
1
‍
f
U
1
(
u
)
d
u
∫
k
∞
‍
f
U
1
(
u
)
d
u
.
With the help of a mathematical software, we can obtain numerical values of this probability for any density function f
U
1 and any value of k.Next, ifT is exponentially distributed, it is well known that, given that there is a signal (i.e., an event) between t and t
+
1, the conditional density of T is uniform on the interval (
t
,
t
+
1
], so that the signal in question occurred on average at time t
+
1
/
2. In the general case, we must compute the conditional density function of U
1 in the interval (
k
,
k
+
1
], given that U
1
>
k.If there is indeed a signal betweent and t
+
1, the forecast of the river flow at time t
+
1 would be
(10)
Z
k
:
=
X
(
t
)
e
-
1
/
c
^
+
E
[
Y
]
1
P
[
U
1
>
k
]
×
∫
k
k
+
1
‍
e
-
(
k
+
1
-
u
)
/
c
^
f
U
1
(
u
)
d
u
.
Therefore, the forecast of the river flow at time t
+
1, given the entire history of the process up to t, is given by
(11)
X
^
(
t
+
1
)
=
X
(
t
)
e
-
1
/
c
^
(
1
-
p
k
)
+
Z
k
p
k
=
X
(
t
)
e
-
1
/
c
^
+
p
k
E
[
Y
]
1
P
[
U
1
>
k
]
×
e
-
(
k
+
1
)
/
c
^
∫
k
k
+
1
‍
e
u
/
c
^
f
U
1
(
u
)
d
u
.


## 2.2. Forecasting the Flow at Timet
+
2

LetU
1
,
k
=
U
1
∣
{
U
1
>
k
}. To forecast the river flow at time t
+
2, given all the information available up to time t, we need the distribution of the random variable S
k
:
=
U
1
,
k
+
U
2, where U
2 is defined in (7).We now assume that the probability that there will be more than two signals in the interval(
t
,
t
+
2
] is negligible (irrespective of the value of k). The density function of S
k is given by the convolution of the density functions of U
1
,
k and U
2:
(12)
f
S
k
(
s
+
k
)
=
∫
k
∞
‍
f
U
1
,
k
(
u
)
f
U
2
(
s
+
k
-
u
)
d
u
.
Since U
2 is a nonnegative random variable, we can write that
(13)
f
S
k
(
s
+
k
)
=
∫
k
s
+
k
‍
f
U
1
,
k
(
u
)
f
U
2
(
s
+
k
-
u
)
d
u
for
s
>
0
.In the case whenT has an exponential distribution with parameter λ, we know that (
S
k
-
k
) has a gamma distribution with parameters α
=
2 and λ. In the other cases, we must generally evaluate the above integral numerically.The value ofX
(
t
+
2
) is given by(14)
X
(
t
+
2
)
=
{
e
-
2
/
c
X
(
t
)
if
U
1
,
k
>
k
+
2
,
e
-
2
/
c
X
(
t
)
+
Y
1
*
e
-
(
k
+
2
-
U
1
,
k
)
/
c
if
U
1
,
k
≤
k
+
2
and
U
2
>
2
-
(
U
1
,
k
-
k
)
,
e
-
2
/
c
X
(
t
)
+
Y
1
*
e
-
(
k
+
2
-
U
1
,
k
)
/
c
if
U
1
,
k
+
U
2
≤
k
+
2
+
Y
2
*
e
-
(
2
-
U
2
)
/
c
(
and
U
3
>
2
-
(
U
1
,
k
-
k
)
-
U
2
)
,where Y
i
* denotes the size of the i
th signal in the interval (
t
,
t
+
2
]. As mentioned above, we assume that
(15)
P
[
U
3
>
2
-
(
U
1
,
k
-
k
)
-
U
2
]
≃
1
.We must first compute(16)
p
1
,
k
:
=
P
[
U
1
,
k
∈
(
k
,
k
+
2
]
]
.
When U
1 has an exponential distribution with parameter λ,
(17)
p
1
,
k
=
P
[
U
1
∈
(
0,2
]
]
=
1
-
e
-
2
λ
.
In general, we have
(18)
p
1
,
k
=
∫
k
k
+
2
‍
f
U
1
(
u
)
d
u
∫
k
∞
‍
f
U
1
(
u
)
d
u
.Next,(19)
p
2
,
k
:
=
P
[
U
1
,
k
+
U
2
∈
(
k
,
k
+
2
]
]
≡
P
[
S
k
∈
(
k
,
k
+
2
]
]
=
∫
0
2
‍
f
S
k
(
s
+
k
)
d
s
=
∫
0
2
‍
[
∫
k
s
+
k
‍
f
U
1
,
k
(
u
)
f
U
2
(
s
+
k
-
u
)
d
u
]
d
s
.Alternatively, using the fact thatU
1
,
k and U
2
(
≥
0
) are independent random variables, we may write that
(20)
p
2
,
k
:
=
P
[
U
1
,
k
+
U
2
∈
(
k
,
k
+
2
]
]
=
∫
k
k
+
2
‍
P
[
U
1
,
k
+
U
2
∈
(
k
,
k
+
2
]
∣
U
1
,
k
=
u
]
f
U
1
,
k
(
u
)
d
u
=
∫
k
k
+
2
‍
P
[
U
2
∈
(
k
-
u
,
k
+
2
-
u
]
]
f
U
1
,
k
(
u
)
d
u
=
∫
k
k
+
2
‍
[
∫
0
k
+
2
-
u
‍
f
U
2
(
v
)
d
v
]
f
U
1
,
k
(
u
)
d
u
.Finally, by independence,(21)
P
[
{
U
1
,
k
≤
k
+
2
}
∩
{
U
2
>
2
-
(
U
1
,
k
-
k
)
}
]
=
∫
k
k
+
2
‍
P
[
U
2
∈
(
k
+
2
-
u
,
∞
)
]
f
U
1
,
k
(
u
)
d
u
=
∫
k
k
+
2
‍
[
∫
k
+
2
-
u
∞
‍
f
U
2
(
v
)
d
v
]
f
U
1
,
k
(
u
)
d
u
=
∫
k
k
+
2
‍
[
1
-
∫
0
k
+
2
-
u
‍
f
U
2
(
v
)
d
v
]
f
U
1
,
k
(
u
)
d
u
=
p
1
,
k
-
p
2
,
k
.The forecast of the river flow at timet
+
2, given the entire history of the process up to time t, is thus given by
(22)
X
^
(
t
+
2
)
=
X
(
t
)
e
-
2
/
c
^
(
1
-
p
1
,
k
)
+
{
X
(
t
)
e
-
2
/
c
^
+
E
[
Y
]
∫
k
k
+
2
‍
e
-
(
k
+
2
-
u
)
/
c
^
f
U
1
,
k
(
u
)
d
u
}
×
(
p
1
,
k
-
p
2
,
k
)
+
{
X
(
t
)
e
-
2
/
c
^
+
E
[
Y
]
∫
k
k
+
2
‍
e
-
(
k
+
2
-
u
)
/
c
^
f
U
1
,
k
(
u
)
d
u
+
E
[
Y
]
∫
k
k
+
2
‍
[
∫
0
k
+
2
-
u
‍
e
-
(
2
-
v
)
/
c
^
f
U
2
(
v
)
d
v
]
kkkkkkkkkkkkk
×
f
U
1
,
k
(
u
)
d
u
[
∫
0
k
+
2
-
u
‍
e
-
(
2
-
v
)
/
c
^
f
U
2
(
v
)
d
v
]
}
p
2
,
k
=
X
(
t
)
e
-
2
/
c
^
+
{
E
[
Y
]
∫
k
k
+
2
‍
e
-
(
k
+
2
-
u
)
/
c
^
f
U
1
,
k
(
u
)
d
u
}
p
1
,
k
+
{
E
[
Y
]
∫
k
k
+
2
‍
[
∫
0
k
+
2
-
u
‍
e
-
(
2
-
v
)
/
c
^
f
U
2
(
v
)
d
v
]
l
l
l
l
l
l
l
l
l
l
×
f
U
1
,
k
(
u
)
d
u
∫
0
k
+
2
-
u
‍
f
U
2
(
v
)
d
v
}
p
2
,
k
=
X
(
t
)
e
-
2
/
c
^
+
{
E
[
Y
]
e
-
(
k
+
2
)
/
c
^
∫
k
k
+
2
‍
e
u
/
c
^
f
U
1
,
k
(
u
)
d
u
}
p
1
,
k
+
{
E
[
Y
]
∫
k
k
+
2
‍
[
e
-
2
/
c
^
∫
0
k
+
2
-
u
‍
e
v
/
c
^
f
U
2
(
v
)
d
v
]
k
k
k
k
k
k
k
k
k
k
k
k
×
f
U
1
,
k
(
u
)
d
u
∫
0
k
+
2
-
u
‍
e
v
/
c
^
f
U
2
(
v
)
d
v
}
p
2
,
k
.


## 3. Forecasting the Flows of the Delaware and Hudson Rivers

To assess the quality of the forecasts obtained with the renewal filtered process, we need to apply the formulas developed in the previous section to real-life data. As in [7], we chose the Delaware and Hudson Rivers, located in the United States, at the Montague (01438500) and the North Creek (01315500) gage stations, respectively. The flow values are available on the Internet at the following address: http://nwis.waterdata.usgs.gov. To calibrate the model, we used the data from October 2008 to September 2009.Remarks. (i) When one is dealing with real-life data, things are not as simple as the mathematical models suggest. Looking at the actual hydrograph of the Delaware River, we observe a number of weak peaks that occur while the flow is decreasing and for which the increasing period lasts only one day before the flow resumes its decline. Similarly, there are sometimes small increases of the flow value occurring just after a minimum was observed that last only one day. In building the data set, there is always a subjective part. We decided to neglect these weak peaks, thus considering only the peaks that appear quite clearly in the hydrograph.(ii) Because the river flow is observed on a daily basis, we must discretize the set of possible values taken by the interarrival timeT. We computed the number k of days elapsed between consecutive peaks in the data set. In the case of the first observed peak, k is the number of days elapsed since the beginning of the calibration period. Then, k is obtained by subtracting the arrival times of the consecutive peaks. Notice that k is greater than or equal to 1.We first consider the Delaware River. For this river, we found that the average valueT
¯ of the interarrival times in the data set is equal to 6.8889 days. Moreover, the standard deviation of the observations is given by s
T
=
3.8977 days. Based on these values, one may at once conclude that T is very unlikely to be exponentially distributed. Indeed, remember that the mean and the standard deviation of an exponential distribution with parameter λ are both equal to 1
/
λ.Next, the value of the point estimate of the (unitless) parameterc (see (6)) is c
^
=
9.2592, and the average value of the magnitude Y of the signals is given by Y
¯
=
984.6121 cubic feet per second.For the distribution of the random variableT, we tested the following models:(i)
an exponential distribution with parameterλ,
(ii)
a lognormal distribution with parametersμ and σ,
(iii)
a gamma distribution with shape parameterα and scale parameter λ,
(iv)
a Weibull distribution with shape parameterλ and scale parameter κ,
(v)
a Rayleigh distribution with parameterσ.Table1 summarizes the results of the chi-square goodness-of-fit tests for each distribution and gives the respective parameter estimates. We find that, as expected, the exponential distribution is rejected, whereas the other distributions are all accepted, according to the P-values of the tests. The empirical cumulative distribution function as well as the distribution functions of the various models considered above is shown in Figure 2. We clearly see that the exponential distribution does not fit the data nearly as well as the other four models.Table 1
Goodness-of-fit tests for the distribution of the random variableT (Delaware River).
Distribution
Exponential
Lognormal
Gamma
Weibull
Rayleigh
Parameters
1/6.8889
1.7574, 0.6187
3.0544, 2.2554
7.7974, 1.8975
5.5817
P-value
0.0344
0.7887
0.7723
0.6927
0.7008Figure 2
Empirical and fitted distribution functions for the random variableT (Delaware River).Now, in the case of the filtered Poisson process (1) with the response function defined in (2), to estimate the river flow at time t
+
1 we compute the expected value of X
(
t
+
1
) given all the past observations. Because of the memoryless property of the exponential distribution, we find that this conditional expectation only depends on the most recent value of the flow, that is, X
(
t
). We obtain (see [7]) that
(23)
E
[
X
(
t
+
1
)
X
]
=
(
t
)
e
-
1
/
c
X
(
t
)
+
E
[
X
(
1
)
]
,
where
(24)
E
[
X
(
1
)
]
=
λ
c
μ
(
1
-
e
-
1
/
c
)
.Similarly, to estimate the river flow at timet
+
2, we only need to calculate
(25)
E
[
X
(
t
+
2
)
∣
X
(
t
)
]
=
E
[
∑
n
:
0
≤
τ
n
≤
t
+
2
‍
Y
n
e
-
(
t
+
2
-
τ
n
)
/
c
∣
X
(
t
)
]
=
e
-
2
/
c
E
[
∑
n
:
0
≤
τ
n
≤
t
‍
Y
n
e
-
(
t
-
τ
n
)
/
c
∣
X
(
t
)
]
+
E
[
∑
n
:
t
<
τ
n
≤
t
+
2
‍
Y
n
e
-
(
t
+
2
-
τ
n
)
/
c
∣
X
(
t
)
]
=
e
-
2
/
c
X
(
t
)
+
E
[
X
(
2
)
]
,
where
(26)
E
[
X
(
2
)
]
=
λ
c
μ
(
1
-
e
-
2
/
c
)
.Moreover, we have that(27)
lim
⁡
t
→
∞
E
[
X
(
t
)
]
=
λ
c
μ
,
(28)
lim
⁡
t
→
∞
Var
⁡
[
X
(
t
)
]
=
λ
c
μ
2
,
(29)
lim
⁡
t
→
∞
ρ
X
(
t
)
,
X
(
t
+
δ
)
=
e
-
δ
/
c
,
where ρ
X
(
t
)
,
X
(
t
+
δ
) denotes the correlation coefficient of X
(
t
) and X
(
t
+
δ
).By making use of these three equations, we can estimate the parameterλ of the Poisson process, the parameter μ of the assumed exponential distribution of the Y
n’s, and the parameter c in the response function. Actually, we do not need point estimates of λ and μ to forecast the flow. From (27), we deduce that we can simply estimate the ratio λ
c
/
μ that appears in both (24) and (26) by the mean value of the observed data when the process is in steady state. Finally, to estimate the parameter c, we first compute the empirical correlation coefficient r
1 between the flow values at times t and t
+
1 and then we set c
^
=
-
1
/
ln
⁡
(
r
1
) (see (29)). With the value of r
1 being close to 1, the previous formula is well defined.Remark. If one is only interested in forecasting the river flow at time t
+
δ, one may estimate c by computing the empirical value of ρ
X
(
t
)
,
X
(
t
+
δ
) and making use of (29).We are now ready to compare the forecasts derived from both the filtered Poisson process and the filtered renewal process with the various distributions considered for the interarrival timeT. We used the formulas given above to forecast the flow of the Delaware River for the 89-day period from February to April 2010.To assess the accuracy of the forecasts, we consider four criteria commonly used in hydrology: themean absolute percentage error (Mape), theNash criterion (Nash), thepeak criterion (Pc), and thecorrelation coefficient (Corr) between the observed and forecasted values.The Nash criterion is based on the mean square error and is used to assess the predictive power of hydrological models. It evaluates the quality of the forecasts from the differences between the expected and observed daily values. It is equal to 1 in case of a perfect fit. For its part, the peak criterion is used to measure the quality of the forecasts during the critical peak period. The closer to 0 it is, the better the forecasts are.Table2 displays the results obtained with the two competing models. Looking at these results, we notice that the Mape, Nash, and Corr criteria yield practically the same values for the filtered renewal and filtered Poisson processes. However, the filtered renewal process did much better than the filtered Poisson process when we consider the peak criterion. This criterion is especially important because accurate forecasts are really needed during the peak period.Table 2
Criteria for comparing the forecasting models (Delaware River).
Filtered renewal process
Filtered Poisson process
Lognormal
Gamma
Weibull
Rayleigh
Exponential
X
(
t
+
1
)
Mape
12.58%
12.49%
12.55%
12.56%
12.76%
Nash
0.8288
0.8290
0.8287
0.8287
0.8356
Pc
0.1674
0.1664
0.1663
0.1669
0.2522
Corr
0.9188
0.9186
0.9180
0.9188
0.9187
X
(
t
+
2
)
Mape
20.36%
20.33%
20.34%
20.34%
21.83%
Nash
0.5447
0.5445
0.5445
0.5445
0.5573
Pc
0.2930
0.2893
0.2869
0.2863
0.3344
Corr
0.7735
0.7732
0.7731
0.7730
0.7734We now turn to the Hudson River. The value of the point estimate of the constantc for this river is c
^
=
9.1172 and the average magnitude of the signals is Y
¯
=
229.2533, which is much smaller than in the case of the Delaware River.Proceeding as above, we performed chi-square goodness-of-fit tests for the various distributions considered. The results are presented in Table3 (see also Figure 3). This time, we see that both the exponential and Rayleigh distributions must be rejected, while the lognormal distribution clearly provides the best fit to the data.Table 3
Goodness-of-fit tests for the distribution of the random variableT (Hudson River).
Distribution
Exponential
Lognormal
Gamma
Weibull
Rayleigh
Parameters
1/8.3488
1.8754, 0.6625
9.2134, 1.3399
2.1785, 3.8324
7.9978
P-value
0.0013
0.6983
0.3203
0.2532
0.000139Figure 3
Empirical and fitted distribution functions for the random variableT (Hudson River).The values of the four criteria used to compare the models are shown in Table4. We discarded the Rayleigh distribution, because it yielded a very bad fit to the data. As in the case of the Delaware River, the values of the Nash and Corr criteria are quite similar, while the important peak criterion and the Mape criterion are much smaller for the filtered renewal process than for the filtered Poisson process.Table 4
Criteria for comparing the forecasting models (Hudson River).
Filtered renewal process
Filtered Poisson process
Lognormal
Gamma
Weibull
Exponential
X
(
t
+
1
)
Mape
12.33%
12.14%
12.33%
15.58%
Nash
0.8240
0.8243
0.8240
0.8360
Pc
0.2875
0.2895
0.2872
0.4818
Corr
0.9146
0.9146
0.9146
0.9146
X
(
t
+
2
)
Mape
20.84%
21.00%
21.15%
32.78%
Nash
0.5841
0.5838
0.5835
0.6222
Pc
0.4558
0.4533
0.4518
0.5575
Corr
0.7908
0.7909
0.7910
0.7910Finally, we present the observed and forecasted values of the flows of the Delaware and Hudson Rivers at timest
+
1 and t
+
2 in Figures 4, 5, 6, and 7, based on the filtered Poisson process (FPP) and the filtered renewal process (FRP) with a lognormal distribution for the interarrival time T. Notice that the filtered renewal process is generally better than the filtered Poisson process at forecasting high flow values, which is our main concern.Figure 4
Observed and forecasted values of the flow of the Delaware River at timet
+
1.Figure 5
Observed and forecasted values of the flow of the Delaware River at timet
+
2.Figure 6
Observed and forecasted values of the flow of the Hudson River at timet
+
1.Figure 7
Observed and forecasted values of the flow of the Hudson River at timet
+
2.


## 4. Conclusion

In this paper, we were able to significantly improve the short-term hydrological forecasts produced by filtered Poisson processes by using a more realistic distribution for the interarrival time of the signals, thus working with filtered renewal processes. As we saw in the previous section, the filtered Poisson process, despite its lack of realism, is able to yield reasonable short-term forecasts. One possible explanation is the fact that the forecasting formulas derived from this model are very easy to implement. Also, contrary to the filtered renewal process, it does not require subjective decisions to estimate its various parameters. However, the filtered Poisson process was not able to produce forecasts as accurate as those derived from the filtered renewal process during the important peak period.In addition to being used to forecast flow values, the models considered in this paper can also give us an estimate of the probability that the river flow will exceed a certain threshold in the next few days. This threshold may be a value that corresponds to a flow above which the risk of flooding is very high.Finally, we could use the response function in (3) to further improve the forecasts. To do so, we would have to be able to estimate the parameter d in that response function and then obtain formulas that generalize the ones in (11) and (22), which is probably not an easy task.



---

*Source: 593243-2014-05-05.xml*