Estimating the index of increase via balancing deterministic and random data
Lingzhi Chen Youri Davydov Nadezhda Gribkova Ričardas Zitikis
math.ST stat.TH 62G05, 62G08, 62G20 (Primary) 62P15, 62P20, 62P25 (Secondary)









Estimating the index of increase via balancing deterministic and random data





Lingzhi Chen^1,[lchen522@uwo.ca],
Youri Davydov^2,[youri.davydov@univ-lille1.fr],
Nadezhda Gribkova^3,[n.gribkova@spbu.ru],
and
Ričardas Zitikis^1,[rzitikis@uwo.ca]




^1School of Mathematical and Statistical Sciences,
Western University, London, Ontario N6A 5B7, Canada



^2Chebyshev Laboratory, St. Petersburg State University, Vasilyevsky Island, St. Petersburg 199178, Russia



^3Faculty of Mathematics and Mechanics, St. Petersburg State University,
 St. Petersburg 199034, Russia







Abstract. We introduce and explore an empirical index of increase that works in both deterministic and random environments, thus allowing to assess monotonicity of functions that are prone to random measurement-errors. We  prove consistency of the index and show how its rate of convergence is influenced by deterministic and random parts of the data. In particular, the obtained results suggest a frequency at which observations should be taken in order to reach any pre-specified level of estimation precision. We illustrate the index using data arising from purely deterministic and error-contaminated functions, which may or may not be monotonic.




Key words and phrases: index of increase, determinism, randomness, measurement errors, smoothing, cross validation.




2010 MSC:  Primary: 62G05, 62G08, 62G20; Secondary: 62P15, 62P20, 62P25.







§ INTRODUCTION




Dynamic processes in populations are often described using functions (e.g., Bebbington et al., 2007, 2011; and references therein). They are observed in the form of data points, usually contaminated by measurement errors. We may think of these points as randomly perturbed true values of underlying functions, whose measurements are taken at certain time instances. The functions, their rates of change, and de/acceleration can be and frequently are  non-monotonic. Nevertheless, it is of interest to assess and even compare the extent of their monotonicity, or lack of it. We refer to Qoyyimi (2015) for a discussion and literature review of various applications.


Several methods for assessing monotonicity have been suggested in the literature (e.g., Davydov and Zitikis, 2005, 2017;  Qoyyimi, and Zitikis, 2014, 2015). In particular, Davydov and Zitikis (2017) show the importance of such assessments in insurance and finance, especially when dealing with weighted insurance calculation principles (Furman and Zitikis, 2008), among which we find such prominent examples as the Esscher (Bühlmann, 1980, 1984), Kamps (1998), and Wang (1995, 1998) premiums. Furthermore, Egozcue et al. (2011) provide problems in economics where the sign of the covariance

    𝐂𝐨𝐯[X,w(X)]

needs to be determined for various classes of function w. One of such examples concerns the slope of indifference curves in two-moment expected utility theory (e.g., Eichner and Wagener, 2009;  Sinn, 1990;  Wong, 2006; and references therein). Another problem concerns decision making (e.g., speculation, normal backwardation, contango, etc.) of competitive companies under price uncertainty (e.g., Holthausen, 1979; Feder et al., 1980; Hey, 1981; Meyer and Robison, 1988; and references therein).


Lehmann (1966) has shown that if the function w is monotonic, then covariance (<ref>) is either positive (when w is increasing) or negative (when w is decreasing). This monotonicity assumption on w, though satisfied in a number of cases of practical interest, excludes a myriad of important cases with more complex risk profiles. For example, when dealing with the aforementioned economics-based problems, the role of w is played by the derivative u' of the underlying utility function, which may not be convex or concave everywhere, as argued and illustrated by, e.g., Friedman and Savage (1948), Markowitz (1952), Kahneman and Tversky (1979), Tversky and Kahneman (1992), among others. Hence, since w might be non-monotonic, how far can this function be from being monotonic, or increasing? Furthermore, since the population risk- or utility-profile cannot be really known, the non-monotonicity of w needs to be assessed from data, and this leads us to the statistical problem of this paper.


In addition, supported by the examples of Anscombe (1973) on potential pitfalls when using the classical correlation coefficient, Chen and Zitikis (2017) argue in favour of using the index of increase, as defined by Davydov and Zitikis (2017), for assessing non-monotonicity of scatterplots. Chen and Zitikis (2017) apply this approach to analyze and compare student performance in subjects such as mathematics, reading and spelling, and illustrate their reasoning on data provided by Thorndike and Thorndike-Christ (2010). One of the methods discussed by Chen and Zitikis (2017) deals with scatterplots representing finite populations, in which case large-sample estimation is not possible. The other method involves large-sample regression techniques (Figure <ref>), in which case Chen and Zitikis (2017) calculate the corresponding indices of increase using a numerical approach, that gives rise to the values denoted by I and reported in the bottom-right corners of the panels of Figure <ref>.

Though important, these methods do not allow direct large-sample  non-monotonicity  quantifications and thus inferences about larger populations. In this paper, therefore, we offer a statistically attractive and computationally efficient procedure for assessing data patterns that arise from non-monotonic patterns contaminated by random measurement errors.

We have organized the rest of the paper as follows. In Section <ref>, we introduce the index and provide basic arguments leading to it. In Section <ref>, we explain why and how the index needs to be adjusted in order to become useful in situations when random measurement errors are present. In Section <ref>, we rigorously establish consistency of the estimator and introduce relevant data-exploratory and cross-validatory techniques.  Since the limiting distribution of the estimator is complex, in Section <ref> we implement a bootstrap-based procedure for determining standard errors and, in turn, for deriving confidence intervals. Section <ref> concludes the paper with a brief summary of our main contributions.





§ THE INDEX OF INCREASE




Davydov and Zitikis (2017) have introduced the index of increase

    I(h_0)=∫_a^b (h_0')_+dλ∫_a^b |h_0'|dλ   (:=∫_a^b (h_0'(t))_+dt ∫_a^b |h_0'(t)|dt  )

for any absolutely continuous (e.g., differentiable) function h_0 on interval [a,b], where (h_0')_+:=max{h_0' , 0}, and “:=” denotes equality by definition. Throughout the paper, we use λ to denote the Lebesgue measure, which helps us to write integrals compactly, as seen from the  ratios above. We shall explain how the index arises later in the current section. Of course, this framework reduces to the unit interval [0,1] by considering the function h(t):=h_0(a+(b-a)t) instead of h_0. Namely, we have

    I(h_0)=∫_0^1 (h')_+dλ∫_0^1 |h'|dλ=:I(h).



To illustrate, in Figure <ref>

we have visualized the following quartet of functions

    h_1(t)=sin(-π 2+3π 2t),
          h_2(t)=cos(-π 2+3π 2t),
    
    
    h_3(t)=sin(π 2t),
          h_4(t)=cos(π 2t),

and we have also calculated their indices of increase. Since h_3 and h_4 are monotonic functions on the interval [0,1], calculating their indices of increase using formula (<ref>) is trivial, but the same task in the case of non-monotonic functions h_1 and h_2 requires some effort. To facilitate such calculations in a speedy fashion, and irrespective of the complexity of functions, we suggest using the numerical approximation

    I_n(h):=∑_i=2^n (h(t_i,n)-h(t_i-1,n))_+∑_i=2^n |h(t_i,n)-h(t_i-1,n)|

with t_i,n=(i-1)/(n-1) for i=1,…, n. Intuitively, I_n(h) is the proportion of the upward movements of the function h with respect to all the movements, upward and downward.


Knowing the convergence rate of I_n(h) to I(h) when n→∞ is important as it allows us to set a frequency n at which the measurements of h(t_i,n) could be taken during the observation period (e.g., unit interval [0,1]) so that any pre-specified estimation precision of I(h) would be achieved. For example, we have used n=10000 to calculate the index values with the four-digit precision reported in Figure <ref>. We refer to Chen and Zitikis (2017) for details on computational precision.

The following proposition, which is a special case of Lemma <ref> below, establishes the convergence rate based on the level of smoothness of the function h.




Let h be a differentiable function defined on the unit interval [0,1], and let its derivative h' be γ-Hölder continuous for some γ∈ (0,1]. Then, when n→∞,  we have

    ∑_j=2^n ℓ (h(t_i,n)-h(t_i-1,n) )
    = ∫_0^1ℓ ( h'  )dλ + O(n^-γ)

for any positively homogeneous and Lipschitz function ℓ (e.g., ℓ(t)=t_+ and ℓ(t)=|t|). Consequently,

    I_n(h)=I(h)+O(n^-γ).





To explain the basic meaning of the index I(h), we start with an un-normalized version of it, which we denote by J(h). Namely, let ℱ denote the set of all absolutely continuous functions f on the interval [0,1] such that f(0)=0. Denote the total variation
of f∈ℱ on the interval [0,1] by ‖ f ‖, that is, ‖ f ‖ =∫_0^1 |f'|dλ. Furthermore, by definition, we have (f')_+=max{f' , 0} and (f')_-=max{ -f', 0}, and we also have the equations f'=(f')_+-(f')_- and  |f'|=(f')_++(f')_-. Finally, we use ℱ^- to denote the set of all the functions f∈ℱ that are non-increasing. All of these are of course well-known fundamental notions of Real Analysis (e.g., Kolmogorov and Fomin, 1970; Dunford and Schwartz,1988; and Natanson, 2016).


For any function h∈ℱ, we define its (un-normalized) index of increase J(h) as the distance between h and the set ℱ^-, that is,

    J(h)
    =inf_f∈ℱ^-‖ h-f‖ .

Obviously, if h is non-increasing, then J(h)=0, and the larger the value of J(h), the farther the function h is from being non-increasing on the interval [0,1]. Determining the index J(h) using its definition (<ref>) is not, however, a straightforward task, and to facilitate it, we next establish a very convenient integral representation of J(h).



The infimum in definition (<ref>) is attained at any function
f_1∈ℱ^- such that f_1'=-(h')_-, and thus

    J(h)=∫_0^1 (h')_+dλ .



A direct proof of this theorem was not provided by Davydov and Zitikis (2017), who refer to a more general and abstract result. Nevertheless, a short and enlightening proof exists, and we present it next.


We start with the note that the bound J(h)≤‖ h - f ‖ holds for every
function f∈ℱ^-, and in particular for the function
f_1 specified in the formulation of the theorem. Hence,

    J(h)
       ≤∫_0^1 |h'-f_1'|dλ
       = ∫_0^1 |h'+(h')_-|dλ
       = ∫_0^1 (h')_+dλ .

It now remains to show the opposite bound. Let T^+ be the set of all
t∈ [0,1] such that h'(t)> 0, and let  T^- be the complement of
the set  T^+, which consists of all those t∈ [0,1] for which
h'(t)≤ 0. Then

    J(h)
       =inf_f∈ℱ^- ( ∫_T^+ |h'-f'|dλ
    +∫_T^- |h'-f'|dλ )
    
       ≥inf_f∈ℱ^-∫_T^+ |h'-f'|dλ
       = inf_f∈ℱ^- (∫_T^+ h'dλ
    + ∫_T^+ |f'|dλ )
    
       = ∫_0^1 (h')_+dλ ,

where the last equation holds when f'(t)=0 for all t∈ T^+, that is, when f'=-(h')_-.  Bounds (<ref>) and (<ref>) establish equation (<ref>), thus finishing the proof of Theorem <ref>.




The index J(h) never exceeds ‖ h‖, and so the normalized version of J(h) is

    I(h):=J(h)/‖ h‖,

which is exactly the index of increase given by equation (<ref>). In summary, the index of increase I(h) is the normalized distance of the function h from the set ℱ^- of all non-increasing functions on the interval [0,1]: we have  I(h)=0 when the function h is non-increasing, and I(h)=1 when the function is non-decreasing. The closer the index I(h) is to 1, the more (we say) the function h is increasing, and the closer it is to 0, the less (we say) the function h is increasing or, equivalently, the more it is decreasing.





§ PRACTICAL ISSUES AND THEIR RESOLUTION



Measurements are usually taken with errors, whose natural model is some distribution (e.g., normal) with mean 0 and finite variance σ^2. In other words, the numerical index I_n(h) turns into the random index of increase

    I_n(h,ε):=∑_i=2^n (Y_i,n-Y_i-1,n)_+∑_i=2^n |Y_i,n-Y_i-1,n| ,

where, for i=1,…, n,

    Y_i,n=h(t_i,n)+ε_i .



Right at the outset, however, serious issues arise. To illustrate them in a speedy and transparent manner, we put aside mathematics such as in Davydov and Zitikis (2004, 2007) and, instead, simulate n=10000 standard normal errors ε_i, thus obtaining four sequences Y_i,n corresponding to the functions of quartet (<ref>).  Then we calculate the corresponding indices of increase using formula (<ref>). All of the obtained values of I_n(h) are virtually equal to 1/2 (see Figure <ref>).

Clearly, there is something amiss.


It is not, however, hard to understand the situation: when all ε_i's are zero, the definition of the integral as the limit of the Riemann sums works as intended, but when the ε_i's are not zero, they accumulate so much when n gets larger that the deterministic part (i.e., the Riemann sum) gets hardly, if at all, visible (compare Figures <ref> and <ref>). In summary, we are facing two extremes:

  
  * If the model is purely deterministic in the sense that there are no measurement errors, which we can understandably argue to be outside the realm of practice, then the more frequently we observe the function h, the more precisely we can estimate its index of increase.
  
  * If, however, there are measurement errors, as they usually are in practice, then the more frequently we observe the function, the less precisely we can estimate   its index of increase, because the accumulated measurement errors obscure the deterministic part.

Neither of the two extremes can be of much interest, or use, for reasons either practical or computational. The purpose of this paper is to offer a way out of this difficulty by showing how to strike a good balance between determinism and randomness inherent in the problem.

We next present an intuitive consideration that will guide our subsequent mathematical considerations, and it will also hint at potential applications of this research. Namely, suppose that the unit interval [0,1] represents an one-day observation period, and let an observation be taken (e.g., by a measuring equipment) every second. Hence, in total, we have n=86400 observations Y_i,n of the (unknown) function h, and they are prone to measurement errors ε_i as in expression (<ref>). For the sake of argument, let ε_i's be i.i.d. standard normal. If we calculate the index I_n(h,ε) based on these data, we already know the problem: I_n(h,ε) tends to 1/2 when n→∞. To diminish the influence of these errors, we average the observed values:

    1 n∑_i=1^n Y_i,n   =1 n∑_i=1^n h(t_i,n)+1 n∑_i=1^n ε_i
    
       d≈∫_0^1 hdλ +1√(n)ε_0,

where d≈ means  `approximately in distribution,' and ε_0 follows the standard normal distribution. However, in the process of averaging out the errors, we have inevitably also averaged the deterministic part and arrived at the mean value ∫_0^1 hdλ of the function h. This value has very little to do with the index I(h), which fundamentally relies on the derivative h'. In short, we have clearly over-averaged the observations Y_i,n: having maximally reduced the influence of measurement errors, we have obscured the function h so much that the estimation of I(h) has become impossible. Clearly, we need to adopt a more tempered approach.

Hence, we group the observations into only M<n groups G_j,n, j=1,… , M, whose cardinalities N:=#(G_j,n) we assume to be the same for all j=1,… , M. It is convenient to re-parametrize these choices using parameter α∈ (0,1), which turns M and N into

    M=⌊ n^α⌋  and   N=⌊ n^1-α⌋ .

This re-parametrization is not artificial. It is, in a way, connected to smoothing histograms and estimating regression functions, and in particular to bandwidth selection in these research areas. We shall elaborate on this topic more in the next section. At the moment, we only note that the aforementioned connection plays a pivotal role in obtaining practically useful and sound estimates of the parameter α.

To gain additional intuition on the grouping parameter α, we come back for a moment to our numerical example with the one-day observation period, which is comprised of n=86400 observations, one per second. Suppose that we decide to average the sixty observations within each minute. Thus, we have N=60 and in this way produce M=1440 new data points, which we denote by Y_j,n. Since NM=n, we have α=1-log(N)/log(n) and thus α =0.6398. If, however, instead of averaging minute-worth data we decide to average, for example, hour-worth data, then we have N=360 (=group cardinality), M=240 (=number of groups), and thus α=0.4822.


Continuing our general discussion, we average the original observations Y_i,n, i=1… , n, falling into each group G_j,n and in this way obtain M group-averages

    Y_j,n:=1 N∑_i∈ G_j,n Y_i,n,     j=1,…, M.

Based on these averages, we modify the earlier introduced index I_n(h,ε) as follows:

    I_n,α(h,ε):=∑_j=2^M (Y_j,n-Y_j-1,n)_+∑_j=2^M |Y_j,n-Y_j-1,n| .

The problem that we now face is to find, if exist, those values of α∈ (0,1) that make the index I_n,α(h,ε) converge to I(h) when n→∞. This is the topic of the next section.






§ CONSISTENCY



The following theorem establishes consistency of the estimator I_n,α(h,ε) and, in particular, specifies the range of possible α values.



Let h be a differentiable function defined on the unit interval [0,1], and let its derivative h' be γ-Hölder continuous for some γ∈ (0,1]. If α∈ (0, 1/3), then  I_n,α(h,ε) is a consistent estimator of I(h), that is, when n→∞, we have

    I_n,α(h,ε)𝐏→I(h).

The rate of convergence is of the order

    O_𝐏(1) n^-min{δ(α ), ρ(α ) }

with  δ(α )=αγ arising from the deterministic part of the problem, that is, associated with the function h, and ρ(α )= (1-3α)/2 arising from the random part, that is, associated with the measurement errors ε_i's.


We next discuss the choice of α from the theoretical and practical perspectives, which do not coincide due to a number of reasons, such as the fact that theory is concerned with asymptotics when n→∞, while practice deals with finite values of n, though possibly very large. Under the (practical) non-asymptotic framework, any value of α∈ (0,1] is, in principle, acceptable because the quantities O_𝐏(1) and n^-min{δ(α ), ρ(α ) } in the specification of convergence rate (<ref>) interact, as both of them depend on h and α.


Under the (theoretical) asymptotic framework, the values α=0 and 1 have to be discarded immediately, as we have already noted. The remaining α's should, as Theorem <ref> tells us, be further restricted to only those below 1/3. Since we wish to chose α that results in the fastest rate of convergence, we maximize the function α↦min{δ(α ), ρ(α ) } and get

    α_max=1 3+2γ.

For example, if the second derivative h”(t) is uniformly bounded on the interval [0,1], which is the case in all our illustrative examples, then γ =1 and thus α_max=1/5.



The grouping and averaging technique that we employ is closely related to smoothing in non-parametric density and regression estimation (e.g., Silverman, 1986; Härdle, 1991; Scott, 2015; and references therein). To elaborate on this connection, we recall that the number of groups is M≈ n^α, whose reciprocal

    b:=1/M≈ n^-α

would play the role of `bandwidth.' In non-parametric density and regression estimation, the optimal bandwidth is of the order O( n^-1/5) when n→∞, which in our case corresponds to α_max=1/5. Hence, α =0 means only one bin/group and thus over-smoothing, whereas α =1 means as many bins/groups as there are observations, and thus under-smoothing. Of course, as we have already noted above, the values α=0 and α=1 are excluded, unless all the measurement errors vanish, in which case smoothing is not necessary and thus α =1 can be used, as we indeed did earlier when dealing with the numerical index I_n(h).


Thinking of the role of γ-Hölder continuity of h' on the problem, it is useful to look at two extreme cases: First, when γ =1, we have α_max=1/5 from formula (<ref>), which corresponds (under weak conditions) to the optimal bandwidth O( n^-1/5) in non-parametric density and regression estimation. Second, when no smoothing is applied, like in the case of the histogram density-estimator, then  (under weak conditions) the optimal bandwidth is of the order O( n^-1/3), which corresponds to α_max=1/3 when γ =0, which essentially means boundedness but no continuity of h'.


Hence, choosing an appropriate value of the grouping parameter α is a delicate task. We next discuss two approaches: The first one is data-exploratory (visual) when we assume that we know the population and want to gain insights into what might happen in practice. The second, practice-oriented approach relies on the idea of cross-validation (e.g., Arlot and Celisse, 2010, Celisse, 2008; and references therein) and is designed to produce estimates of α based purely on data.




 §.§ Data exploratory (visual) choice of α



To gain intuition on how to estimate the grouping parameter α from data, we start out with the functions in quartet (<ref>), which we view as populations, and then we contaminate their observations with i.i.d. errors ε_i ∼𝒩(0,1) according to formula (<ref>).


We have visualized the values of the estimator I_n,α(h,ε) with respect to n and α in Figure <ref>,
	
where the hyperplane in each panel is at the height of the corresponding actual index of increase I(h). For each panel, we visually choose a value of α which is in the intersection of the curved surface with the hyperplane, because in this case the index I_n,α(h,ε) is close to the actual index I(h).

Even though the chosen parameter α value, which we denote by α_vi, may not be optimal due to roughness of the surface, it nevertheless offers a sound choice, as we see from Figure <ref>
	
where we depict the convergence of I_n,α(h,ε) to I(h) when n grows. In each panel, the horizontal red `reference' line is at the height of the actual index value.


Note that in panel (a) of Figure <ref>, the visually obtained α_vi=0.35 is slightly larger than 1/3, but we have to say that we had decided on this value (as a good estimate) before we knew the result of Theorem <ref>, and thus before we knew the (theoretical) restriction α<1/3. Nevertheless, we have decided to leave the value α_vi=0.35 as it is, without tempering with our initial guess in any way. As we shall see in next Section <ref>, however, the purely data-driven and based on cross-validation α value is α_cv=0.28, which is within the range (0,1/3) of theoretically acceptable α values.




 §.§ Choosing α based on cross validation




As we have already elucidated, equation (<ref>) connects our present problem with nonparametric regression-function estimation. In the latter area, researchers usually choose the optimal bandwidth as the point at which cross-validation scores become minimal (e.g., Arlot and Celisse, 2010, Celisse, 2008; and references therein). We adopt this viewpoint as well. Namely, given a scatterplot, say (t_i,n,Y_i,n), we cross validate it (computational details and R packages will be described in a moment). Then we find the minimizing value b=b_cv and finally, according to equation (<ref>),  arrive at the `optimal' α_cv via the equation

    α_cv=log(1/b_cv)/log(n).

In Figure <ref>,
	
we see some differences between the values of α_vi and α_cv. Nevertheless, we should not prejudge the situation in any way because in practice, when no hyperplanes can be produced due to unknown values of I(h), only the values of α_cv can be extracted from data. Note, however, that the four values of α_cv reported in the panels of Figure <ref> are in compliance with the condition of Theorem <ref> stipulating that α's must be in the range (0,1/3) in order to have (asymptotic) consistency.


To explore how the grouped estimator I_n,α(h,ε) based on α_cv's actually performs, we have produced Figure <ref>.
	
Naturally, since the respective visual α_vi's and cross-validatory  α_cv's do not coincide, the corresponding values of I_n,α(h,ε) are also different. Which of them are better from the statistical point of view will become clearer only in Section <ref>, where bootstrap-based standard errors and confidence intervals are derived.


We next present a detailed implementation procedure for finding cross-validatory estimates α_cv of the grouping parameter α. Naturally, the help of the R computing language (R Core Team, 2013) becomes indispensable, and we have used a number of R packages to accomplish the task. We also wish to acknowledge the packages  (Wickham, 2009) and  (Sievert et al., 2017) that we have used extensively in this paper to draw two-dimentional plots and interactive surface plots; the latter plots have been pivotal in extracting the values α_vi  visually.


Hence, from the purely practical computational perspective, we now utilize bandwidth selection techniques of kernel-based regression-function estimation in order to get estimates of the grouping parameter. First, for the sake of programming efficiency, we restrict b's to the interval (0.01, 0.99), and we evenly split the latter interval into bins of width (0.99-0.01)/29 ≈ 0.0338, all of which can of course be refined in order to achieve, if desired, smaller computational errors. Hence, from now on, we have thirty equidistant b's, which are b_i≈ 0.01+(i-1)0.0338 for i=1,… , 30. Next we use the common cross-validation method called repeated k-fold cross validation, and we set k=5 for our purpose. The following main steps are:

    
  * 
    For each function h under consideration, we generate n=10000 data points based on equation (<ref>).
	
  * 
We randomly split the given n points into k folds, denoted by D_1,… , D_k,  of roughly equal sizes.
	
  * 
For each value b_i, we use D_1 as the validation set and let other D's be training sets, which we use to fit a kernel regression model. Specifically, we use the function  from the R package , with the parameter  set to , which means that we use the normal kernel. Then we use the validation set D_1 to get the predicted values and calculate one prediction error, defined as the mean-square error and denoted by E_1. We repeat this step until we use up all the folds as our validation sets. Hence, we obtain k prediction errors E_1,… , E_k. Finally, we average these k prediction errors and denote this average by  E_b_i,1.
	
  * 
We repeat Step <ref> for all b_i's, thus arriving at one estimated prediction error for each b_i. Hence, in total, we have E_b_1,1, …, E_b_30,1.
	
  * 
We repeat Steps <ref>–<ref> fifty times, for every b_i, and then take the averages of the corresponding fifty estimated prediction errors. This gives us fifty final estimates, which we denote by E_b_i. For example, for b_1, the final estimate E_b_1 is the average of E_b_1,1, …, E_b_1,50. In summary, after this step, we have E_b_1, … , E_b_30 of the final estimates of the prediction error.
	
  * 
We draw the plot of the b_i's versus the corresponding estimated prediction errors. The b_i that gives the minimal prediction error is denoted by b_cv. Finally, we use equation (<ref>) to get α_cv.








 §.§ Proof of Theorem <ref>



The following lemma, whose special case is Proposition <ref> formulated earlier, plays a pivotal role when proving Theorem <ref>.



Let h be differentiable, and let its derivative h' be γ-Hölder continuous for some γ∈ (0,1]. Furthermore, let ℓ be any positively homogeneous and Lipschitz function. Then there is a constant c<∞ such that, for any set of points s_1:=0<s_2<⋯ <s_M≤ 1,

    ∑_j=2^M ℓ (h(s_j)-h(s_j-1) )
    = ∫_0^s_Mℓ ( h'  )dλ
    + θ c ∑_j=2^M  | s_j-s_j-1 |^1+γ

where θ is such that |θ|≤ 1.




Since ℓ is Lipschitz and h' is γ-Hölder continuous, we have

    ∫_0^s_Mℓ ( h'  )dλ   = ∑_j=2^M ∫_s_j-1^s_jℓ ( h'(s)  )
    -ℓ ( h'(s_j)  )ds
    + ∑_j=2^M  (s_j-s_j-1 ) ℓ ( h'(s_j)  )
    
       = θ c ∑_j=2^M ∫_s_j-1^s_j | s- s_j |^γds
    + ∑_j=2^M  (s_j-s_j-1 ) ℓ ( h'(s_j)  )
    
       = θ c ∑_j=2^M  | s_j -s_j-1|^1+γ
    + ∑_j=2^M  (s_j-s_j-1 ) ℓ ( h'(s_j)  ),

where the values of c<∞ and |θ|≤ 1 might have changed from line to line. Next, we explore the right-most sum of equation (<ref>), to which we add and subtract the right-hand side of equation (<ref>). Then we use the mean-value theorem with some ξ_j∈ [s_j-1,s_j] and arrive at the equations

    ∑_j=2^M  (s_j-s_j-1 ) ℓ ( h'(s_j)  )
       = ∑_j=2^M ℓ (h(s_j)-h(s_j-1) )
    +∑_j=2^M (s_j-s_j-1)  ( ℓ ( h'(s_j)  )
    -ℓ ( h'(ξ_j) )  )
    
       = ∑_j=2^M ℓ (h(s_j)-h(s_j-1) )
    +θ c ∑_j=2^M  | s_j -s_j-1|^1+γ,

where the last equation holds because ℓ is positively homogeneous and Lipschitz, and h' is γ-Hölder continuous. Equations (<ref>) and (<ref>) imply equation (<ref>) and finish the proof of Lemma <ref>.





We start with the equations

    Y_j,n   =1 N∑_i∈ G_j,n  h(t_i,n)
    +1 N∑_i∈ G_j,nε_i
    
       =1 N∑_i∈ G_j,n  h(t_i,n)
    +ε^*_j,n ,

where

    ε^*_j,n=1 N∑_i∈ G_j,nε_i.

We next tackle the deterministic sum on the right-hand side of equation (<ref>), and start with the equation

    1 N∑_i∈ G_j,n  h(t_i,n)
    = n-1 N∑_i=1^N  h (  i-1 n-1+ (j-1)N n-1 )1 n-1

because G_j,n=(j-1)N+{1, … , N} for all j=1,… , M. Consequently,

    1 N∑_i∈ G_j,n  h(t_i,n)
       = n-1 N ( ∑_i=1^N  h (  i-1 n-1+ (j-1)N n-1 )1 n-1 -∫_(j-1)N/(n-1)^jN/(n-1) hdλ )
    
            + n-1 N∫_(j-1)N/(n-1)^jN/(n-1) hdλ
       = n-1 N∫_(j-1)N/(n-1)^jN/(n-1) hdλ  + O(n^-1 ),

where we used the fact that h is Lipschitz. By the mean-value theorem, there is t_j,n^* between (j-1)N/(n-1) and jN/(n-1) such that the right-hand side of equation (<ref>) is equal to h(t_j,n^*)+ O(n^-1).
Consequently, with the notation

    Y^*_j,n:=h(t_j,n^*)+ε^*_j,n,

we have Y_j,n=Y^*_j,n+O(n^-1) and thus the increments  Y_j,n-Y_j-1,n are equal to Y^*_j,n-Y^*_j-1,n+O(n^-1).  This gives us the equations

    ∑_j=2^M ℓ (Y_j,n-Y_j-1,n )
       = ∑_j=2^M ℓ (Y^*_j,n-Y^*_j-1,n ) +O(n^-(1-α))
    
       = ∑_j=2^M ℓ (h(t_j,n^*)-h(t_j-1,n^*) )
    + O ( ∑_j=2^M |ε^*_j,n-ε^*_j-1,n|  )
    +O(n^-(1-α))

because |ℓ(t)-ℓ(s)|≤ |t-s| for all real t and s.
The random variables ε^*_j,n, j=1,…, M, are independent and identically distributed with the means 0 and variances σ^2/N. Hence,

    𝐄 ( ∑_j=2^M |ε^*_j,n-ε^*_j-1,n|  )
       ≤ cM max_j√(𝐄 ( (ε^*_j,n)^2  ) )
       ≤ cM max_j√(σ^2/N )
       =O( n^-(1-3α)/2 ),

which implies

    ∑_j=2^M |ε^*_j,n-ε^*_j-1,n| =O_𝐏(n^-(1-3α)/2 ).

The right-hand side of equation (<ref>) converges to 0 because α∈ (0,1/3). In view of equations (<ref>) and (<ref>), we have

    ∑_j=2^M ℓ (Y_j,n-Y_j-1,n )
    = ∑_j=2^M ℓ (h(t_j,n^*)-h(t_j-1,n^*) )
    + O_𝐏(n^-(1-3α)/2 )+O(n^-(1-α)).

Furthermore, by Lemma <ref> we have

    ∑_j=2^M ℓ (h(t_j,n^*)-h(t_j-1,n^*) )
    = ∫_0^1 ℓ ( h'  )dλ  +O( n^-αγ ).

Combining equations (<ref>) and (<ref>), and using β to denote min{αγ, (1-3α)/2}, we have

    I_n,α(h,ε)
       =∫_0^1 (h')_+dλ +O_𝐏( n^-β ) ∫_0^1 |h'|dλ +O_𝐏( n^-β ) 
       𝐏→∫_0^1 (h')_+dλ∫_0^1 |h'|dλ=I(h) .

The rate of convergence is of the order O_𝐏(n^-β). Theorem <ref> is proved.






§ BOOTSTRAP-BASED CONFIDENCE INTERVALS




To construct confidence intervals for I(h) based on the estimator I_n,α(h,ε), we need to determine standard errors, which turns out to be a very complex task from the viewpoint of asymptotic theory. Hence, we employ bootstrap (e.g., Hall, 1992; Efron and Tibshirani, 1993; Shao and Tu, 1995; Davison and Hinkley, 1997; and references therein). The re-sampling size m is quite often chosen to be equal to the actual sample size n, but in our case, we find it better to re-sample fewer than n observations (i.e., m < n) and thus follow specialized to this topic literature by Bickel et al. (1997), Bickel and Sakov (2008), Gribkova and  Helmers (2007, 2011); see also references therein. Specifically, the steps that we take are:

	
  * For a given function h, we generate n=10000 values y_1, ...., y_n according to the model Y_i = h(t_i,n)+ ε_i, where ε_i are i.i.d. standard normal.
    
  * We re-sample 1000 times and in this way obtain 1000 sub-samples of size m, which we choose to be m ≈ 2√(n) according to a rule of thumb (DasGupta, 2008, p. 478).
	
  * We use formula (<ref>) to calculate the grouped index of increase, thus obtaining 1000 values of it; one value for each sub-sample. We denote the empirical distribution of the obtained values by F^*.
	
  * With Q^* denoting the (generalized) inverse of F^*, the 95% quantile-based confidence interval is (q_2.5%, q_97.5%), where q_2.5%= Q^*(0.025) and q_97.5%= Q^*(0.975).



To illustrate, we introduce a second quartet of functions of this paper, namely:

    h_5(t)=(t-1)^2+ sin(6t),
          h_6(t)=(t-0.25)^2+sin(0.25t),
    
    
    h_7(t)=t^3-5.6t^2+6t,
          h_8(t)=sin(2π t).

We have visualized the functions in Figure <ref>.

As expected, our preliminary analysis has shown that the un-groped estimators converge to 0.5 in all the four cases, but the grouped estimator I_n,α(h,ε) does converge under appropriate choices of the grouping (or smoothing) parameter α values. Next are summaries of our findings using two approaches: the first one is data exploratory (visual) and the second one is based on cross validation.



 §.§ Data exploratory (visual) choice of α


Based on the crossings of surfaces and hyperplanes depicted in Figure <ref>,
	
we choose appropriate α values, denoted by α_vi, for the functions of quartet (<ref>). To check the performance of these values, we draw convergence graphs in Figure <ref>.
	
Next, we use formula (<ref>) to calculate point estimates of the actual index of increase for each of the functions in quartet (<ref>), whose values appear in Table <ref>.

Finally, we use bootstrap to get standard errors and confidence intervals, all of which are also reported in Table <ref>.

Reflecting upon the findings in Table <ref>, we see that the values of α_vi corresponding to the functions h_5 and h_8 are outside the range (0,1/3) specified by the consistency result of Theorem <ref>, but this of course does not invalidate anything – we are simply working with finite sample sizes n. Naturally, we are now eager to compare all the findings reported in Table <ref> with the corresponding ones obtained by cross validation, which is our next topic.




 §.§ Choosing α based on cross validation


We now use the cross-validation technique to get estimates α_cv of the grouping parameter α for all the functions of quartet (<ref>). In Figure <ref>,
	
we visualize the cross-validation scores, specify their minima b_cv, and also report the grouping parameters α_cv derived via the equation α_cv=log(1/b_cv)/log(n). Based on these α_cv values, we explore the performance of I_n,α(h,ε) using the convergence graphs depicted in Figure <ref>.

The values of point estimates, standard errors, and confidence intervals are reported in Table <ref>.


Note that the first three values of α_cv reported in Table <ref> are inside the range (0,1/3) specified by the consistency result of Theorem <ref>, whereas α_cv=0.3396 corresponding to h_8 is just slightly outside the range. Note also that the values of α_cv corresponding to the functions h_5 and h_8 are considerably smaller than the corresponding α_vi's reported in Table <ref>.


The confidence intervals reported in Tables <ref> and  <ref> comfortably cover the actual values of I(h), and the widths of these confidence intervals, denoted by width_vi and width_cv respectively, are comparable for the functions h_5, h_6 and h_7. The width_cv of the cv-based confidence interval for the function h_8 is, however, considerably wider than the corresponding width_vi reported in Table <ref>. In summary, the relative differences width_cv/width_vi-1 for the functions h_5, h_6, h_7 and h_8 are -0.0611, 0.1555, -0.0425 and 1.1541, respectively. We finish the discussion by recalling Wasserman's (2005) advice: “Do not assume that, if the estimator [...] is wiggly, then cross-validation has let you down. The eye is not a good judge of risk” (Remark 20.18, page 317).




§ SUMMARY AND CONCLUDING NOTES




Davydov and Zitikis (2017) introduced an index of increase when populations are modelled with continuous functions. Chen and Zitikis (2017) explored a modification of the index when populations are discrete and presented in the form of scatterplots, and they also explored the situation when scatterplots are viewed as data sets, in which case they fitted (non-monotonic) regression functions and subsequently applied the technique by Davydov and Zitikis (2017) to assess monotonicity of the fitted functions.


In the present paper we have extended the aforementioned technique to the case when it is not desirable, or appropriate, to view scatterplots as populations, or to use regression methods to fit curves to scatterplots. The herein proposed technique is based on grouping and averaging data, and then calculating the index of increase. Since the grouping parameter depends on both deterministic and random features of the underlying problem, we have suggested a way for grouping data so that the resulting estimator of the index of increase would be consistent. Based on this estimator, we have then suggested a construction of bootstrap-based confidence intervals for the index of increase.


The derived theoretical results have been made accessible to practitioners by detailed descriptions and analyses of various computational aspects inherent in our proposed solution of the problem.





§ ACKNOWLEDGEMENTS


We are indebted to the anonymous reviewers for suggestions, insightful comments, and constructive criticism that guided our work on the revision. The research has been supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada.


99


a1973
Anscombe, F.J. (1973).
Graphs in statistical analysis.
American Statistician, 27, 17–21.



ac2010
Arlot, S. and Celisse, A. (2010). A survey of cross-validation procedures for model selection. Statistics Surveys, 4, 40–79.


blz2007
Bebbington, M., Lai, C.D. and Zitikis, R. (2011).
Modeling human mortality using mixtures of bathtub shaped failure distributions.
Journal of Theoretical Biology, 245, 528–538.

blz2011
Bebbington, M., Lai, C.D. and Zitikis, R. (2011).
Modelling deceleration in senescent mortality.
Mathematical Population Studies, 18, 18–37.

bgz1997
Bickel, P.J., Götze, F. and van Zwet, W.R. (1997).
Resampling fewer than  n observations: gains, losses, and remedies for losses.
Statistica Sinica, 7, 1–31.

bs2008
Bickel, P.J. and Sakov, A. (2008).
On the choice of m in the m out of n bootstrap and confidence bounds for extrema.
Statistica Sinica, 18, 967–985.



b1980
Bühlmann, H. (1980).
An economic premium principle.
ASTIN Bulletin, 11, 52–60.



b1984
Bühlmann, H. (1984).
The general economic premium principle.
ASTIN Bulletin, 14, 13–21.


ac2010
Celisse, A. (2008).
Model Selection via Cross-Validation in Density Estimation, Regression, and Change-Points Detection. Université Paris Sud – Paris XI, Paris.
HAL Id: tel-00346320. 


cz2017
Chen, L. and Zitikis, R. (2017).
Measuring and comparing student performance: a new technique for assessing  directional associations. Education Sciences, 7, 1–27.

da2008
DasGupta, A. (2008).
Asymptotic Theory of Statistics and Probability.
Springer, New York.



dh1997
Davison, A.C. and Hinkley, D. V. (1997).
Bootstrap Methods and their Application. Cambridge University Press, Cambridge, UK.


dz2004
Davydov, Y. and Zitikis, R. (2004).
The influence of deterministic noise on empirical measures
generated by stationary processes.
Proceedings of the American Mathematical Society, 132, 1203–1210.


dz2005
Davydov, Y. and Zitikis, R. (2005).
An index of monotonicity and its estimation: a step
beyond econometric applications of the Gini index.
Metron,
63 (special issue in memory of Corrado Gini), 351–372.

dz2007
Davydov, Y. and Zitikis, R. (2007).
Deterministic noises that can be statistically
distinguished from the random ones.
Statistical Inference for Stochastic Processes, 10, 165–179.


dz2017
Davydov, Y. and Zitikis, R. (2017).
Quantifying non-monotonicity of functions and the lack of positivity in signed measures.
Modern Stochastics: Theory and Applications, 4, 219–231.



ds1988
Dunford, N. and Schwartz, J.T. (1988).
Linear Operators, Part 1: General Theory. Wiley, New York.


et1993
Efron, B. and Tibshirani, R.J. (1993). An Introduction to the Bootstrap.  Chapman and Hall/CRC, Boca Raton, FL.


efwz2011
Egozcue, M., Fuentes García, L., Wong, W.K. and Zitikis, R.  (2011).
The covariance sign of transformed random variables with applications to economics and finance.
IMA Journal of Management Mathematics, 22, 291–300.



ew2009
Eichner, T. and Wagener, A. (2009).
Multiple risks and mean-variance preferences.
Operations Research, 57, 1142–1154.




fjs1980
Feder, G., Just, R.E. and Schmitz, A.  (1980).
Futures markets and the theory of the firm under price uncertainty.
Quarterly Journal of Economics, 94, 317–328.




fs1948
Friedman, M. and Savage, L.J. (1948).
The utility analysis of choices involving risk.
Journal of Political Economy, 56, 279–304.



fz2008
Furman, E. and Zitikis, R. (2008).
Weighted premium calculation principles.
Insurance: Mathematics and Economics, 42, 459–465.


gh2007
Gribkova, N.V. and  Helmers, R. (2007).
On the Edgeworth expansion and the M out of  N bootstrap accuracy  for a Studentized trimmed mean.
Mathematical Methods of Statistics, 16, 142–176.

gh2011
Gribkova, N.V. and Helmers, R. (2011).
On the consistency of the M ≪ N bootstrap approximation for a trimmed mean.
Theory of Probability and Its Applications, 55, 42–53.



h1992
Hall, P. (1992).
The Bootstrap and Edgeworth Expansion. Springer, New York.


h1991
Härdle, W. (1991).
Smoothing Techniques, with Implementation in S.
Springer, New York.



h1981
Hey, J.D. (1981).
Hedging and the competitive labor-managed firm under price
uncertainty hedging and the competitive labor-managed firm under price uncertainty.
American Economic Review, 71, 753–757.



h1979
Holthausen, D.M. (1979).
Hedging and the competitive firm under price uncertainty hedging and
the competitive firm under price uncertainty.
American Economic Review, 69, 989–995.



k1998
Kamps, U. (1998).
On a class of premium principles including the Esscher principle.
Scandinavian Actuarial Journal, 1998, 75–80.



kt1979
Kahneman, D. and Tversky, A. (1979). Prospect theory of
decisions under risk. Econometrica, 47, 263–291.


kf1970
Kolmogorov, A.N. and Fomin, S.V. (1970). Introductory Real Analysis.
Dover, New York.



l1966
Lehmann, E.L. (1966). Some concepts of dependence.
Annals of Mathematical Statistics, 37, 1137–1153.



m1952
Markowitz, H. (1952).
The utility of wealth.
Journal of Political Economy, 60, 151–156.



mr1988
Meyer, J. and Robison, L.J. (1988).
Hedging under output price randomness hedging under output price randomness.
American Journal of Agricultural Economics, 70, 268–272.


n2016
Natanson, I.P. (2016).
Theory of Functions of a Real Variable. Dover, New York.


q2015
Qoyyimi, D.T. (2015). A Novel Method for Assessing Co-monotonicity: an Interplay between Mathematics and Statistics with Applications. Electronic Thesis and Dissertation Repository Nr. 3322. 

qz2014
Qoyyimi, D.T. and Zitikis, R. (2014).
Measuring the lack of monotonicity in functions.
Mathematical Scientist, 39, 107–117.

qz2015
Qoyyimi, D.T. and Zitikis, R. (2015).
Measuring association via lack of co-monotonicity: the LOC index and a problem of educational assessment.
Dependence Modeling, 3, 83–97.

r2013
R Core Team (2013). R: A Language and Environment for Statistical Computing.
R Foundation for Statistical Computing, Vienna.


cctskmp2017
Sievert, C., Parmer, C., Hocking, T., Chamberlain, S., Ram, K.,
Corvellec, M. and Despouy, P. (2017).
plotly: Create Interactive Web Graphics via 'plotly.js'. 

s1986
Silverman, B.W. (1986). Density Estimation for Statistics and Data Analysis. Chapman and Hall/CRC, London.

s2015
Scott, D.W. (2015).
Multivariate Density Estimation: Theory, Practice, and Visualization. (Second edition.)
Wiley, New York.

st1995
Shao, J. and Tu, D. (1995).
The Jackknife and Bootstrap.  Springer, New York.


s1990
Sinn, H.-W. (1990).
Expected utility, μ-σ preferences, and linear
distribution classes: a further result.
Journal of Risk and Uncertainty, 3, 277–281.



tt2010
Thorndike, R.M. and Thorndike-Christ, T. (2010).
Measurement and Evaluation in Psychology and Education (8th edition). Prentice Hall, Boston, MA.


tk1992
Tversky, A. and Kahneman, D. (1992). Advances in prospect
theory: Cumulative representation of uncertainty. Journal of
Risk and Uncertainty, 5, 297–323.


w1995
Wang, S. (1995).
Insurance pricing and increased limits ratemaking by proportional
hazards transforms.
Insurance: Mathematics and Economics, 17, 43–54.



w1998
Wang, S. (1998).
An actuarial index of the right-tail risk.
North American Actuarial Journal, 2, 88–101.

w2005
Wasserman, L. (2005).
All of Statistics: a Concise Course in Statistical Inference. (Second printing). Springer, New York.

w2009
Wickham, H. (2006). ggplot2: Elegant Graphics for Data Analysis;
Springer, New York.



w2006
Wong, W.K. (2006).
Stochastic dominance theory for location-scale family.
Journal of Applied Mathematics and Decision Sciences, 2006, 1–10.







