Approximating a Complementary Cumulative Distribution Function via a piece-wise function

I hope this is not too much to read, but I tried to give you a specific overview over my problem.

I am currently trying to model the German electricity market, with a special focus on balancing power.

Considering this setting I have a data set consisting of hourly observations of balancing power in Megawatts for the year of 2013. As the application of this data set would not represent any uncertainty in the data I wanted to use a little trick to "assume" uncertainty. While this is really only a basic approximation my model has to have some sort of "uncertainty", else there will be a big mistake. In order to do this, I divided my data set into five blocks (too many more would be to much for my calculation) of varying Megawatt-steps calculating the "probability" of a specific value to be in one of these blocks using the Complementary cumulative distribution function.

What I mean by that is that I took the aforementioned function and tried to best approximate it in a (5-)piece-wise fashion. Unfortunately I had to do this in my example by a "sense of proportion". I hope the picture makes it a little more clear (On the x-axis the values are in Megawatts).

EDIT: To be clear, I am not trying to best-possibly approximate the whole-function in a sense of a non parametric boundary estimation as user603 kindly described in his answer. (Or maybe this can be done using non parametric boundary estimation...?) I just want to fit the five (or whatever number) blocks best-possible to the function in a fashion that, e.g. the Sum of squared Errors is at a minimum. Also: i would love the knot location to be estimated.

The following code example should further explain what I do:

# This function divides the input data into 5 different blocks.
# v, w, x and y are all Megawatt values dividing the data set
getCall_block <- function(v, w, x, y, data){

n <- nrow(data)
#' parameters for function
data_v <- subset(data, data[,2] <= v, select = c(colnames(data)[1],
colnames(data)[2]))
data_w <- subset(data, data[,2] > v & data[,2] <= w, select = c(colnames(data)[1],
colnames(data)[2]))
data_x <- subset(data, data[,2] > w & data[,2] <= x, select = c(colnames(data)[1],
colnames(data)[2]))
data_y <- subset(data, data[,2] > x & data[,2] <= y, select = c(colnames(data)[1],
colnames(data)[2]))
data_r <- subset(data, data[,2] > y, select = c(colnames(data)[1],
colnames(data)[2]))
#'
df_v <- getAverage(data, data_v, 0, v)
df_w <- getAverage(data, data_w, v, w)
df_x <- getAverage(data, data_x, w, x)
df_y <- getAverage(data, data_y, x, y)
df_r <- getAverage(data, data_r, y, max(data[,2]))

return(rbind(df_v, df_w, df_x, df_y, df_r))
}


And the getAverage function...

#This function is called in the above mentioned function.
#It calculates the mean average of the points of the CCDF lying in one block.
#v and w are the input parameters in Megawatts.
#They can be understood as starting point and end point of a given block.

getAverage <- function(data_c, data, v, w){
#' Imposing greater equal on v if v > 0
if(w > v){
if (v > 0){
v <- v + 1
} else {
v <- 0
}
#' Initializing div, df, mean, n, sum
div <- 0
mean <- 0
n  <- nrow(data)
sum <- 0
df <- data.frame(ID = data[,1], call_freq = numeric(n))
#' Creating inverse ECDF based on dataset
F1 <- ecdf(data_c[,2])
#' Building sum over every data point in the interval [v;w]/(v;w]
for (i in seq(v, w, 1)){
sum <- sum + (1 - F1(i))
div <- div + 1
}
mean <- sum / div
df$call_freq = mapply(function(x) return (x), mean) return (merge(data, df, by = "ID")) } else { cat("Check input parameters. v is greater than w.") } }  My questions are the following: I) Does this way of approximating a CCDF sound okay for you guys? What did I overlook or completely mess up? II) Is there a way to automatically approximate the block sizes, in order to not rely on a sense of proportion. E.g. minimizing the sum of squared errors between the piece-wise function and the original function. Maybe this could be done by using a clustering algorithm like k-means? But then how to cluster this CCDF...? III) Is there a far more easier way to approximate any function in a piece-wise fashion? Thanks in advance! 1 Answer First, It's not really clear to me what you are trying to do. Your question is indeed too long winded and I think this is preventing other users on our site from offering you more answers to choose from. Anyway, I will take a stab at it, hopefully it'll help restructure you question. In the convention, tools and approaches I will describe I will assume that the function you try to approximate is monotone increasing. In your case this means working with$1-\hat{F}$(1 minus the empirical cdf) but that shouldn't be a problem. Now, to answer your questions: I) Does this way of approximating a CCDF sound okay for you guys? What did I overlook or completely mess up? No: you should add more constraints to your optimization problem. For example, you know --through elementary subject knowledge-- that the CDF (CCDF) is monotonous increasing (decreasing). You should explicitely constraint the approximations to also be. Including this constraint will improve the finite sample performance of the estimates and will lessen the dependence of the result on the choice of the nkot points. As to how to do so, I explain in the answer to your question iii) below. Also, you probably want a lower (upper) bound approximation to the (CCDF) CDF for all$x<x_M$where$x_M\$ is your estimate of the largest possible observation. You should also add this as a constraint in the fitting procedure.

II) Is there a way to automatically approximate the block sizes, in order to not rely on a sense of proportion. E.g. minimizing the sum of squared errors between the piece-wise function and the original function. Maybe this could be done by using a clustering algorithm like k-means? But then how to cluster this CCDF...?

If you consider the knot location as fixed, yes: this is a convex optimization problem. If you consider the knot locations as a parameter, then the resulting problem is not convex but can be approximated by a convex optimization problem. In principle, this could be done by adapting the approach I suggest below --in the answer to your question iii)-- in much the same way as is done to transform B splines into penalized spline. This is because the problem you try to solve can be framed as a form of (constrained!) smoothing spline, whose objective function can then be further penalized, to force the algorithm to select the optimal knot points, in much the same way as it done in the case of P-splines. Again, see the answer to iii) below for details.

III) Is there a far more easier way to approximate any function in a piece-wise fashion?

A branch of statistics called non parametric boundary estimation deals with this and related problems (at least to the extend that I understand your question). The idea is to fit piece wise function to your data under the constraint that:

• these piece wise function joint at pre-specified knot points,
• The fitted function is always below the real function evaluated at the knots,
• have a positive derivative (are increasing)
• (possibly) have negative second derivatives (are concave).

There are different methods depending on whether you want the piece-wise function to be linear, quadratic or cubic functions. For the last two, you can also decide whether you want to impose concavity or not (by imposing the fourth constraints as well or only the first three). Obviously, the numerical complexity of the approach depends on how many of these constraints you wish to impose but in any case it involves solving a convex optimization problem (linear program whose size depends linearly on the number of knots for the first three constraints, cone programming problem of size depending (linearly?) on the number of knot for the last constraint).

You will find many links and R implementation at the benchmarking and npbr packages. Here is a simple example comparing DEA, FDH and LFDH (three approach in increasing order of complexity). DEA appears very close to what you are doing already, with the notable difference that it uses general trapezoid instead of box (though it can also use boxes), that in the package I link to, the code is implemented in C++ (so it will probably be a notable speed up over yours) and that the constraints I mentioned above are taken into account. Here is a simple example to illustrate (the original CDF is in black):

library(npbr)

b0<-rexp(1000)
b1<-ecdf(b0)
x<-seq(min(b0),max(b0),l=100)
b3<-b1(x)
plot(x,b3,type="l")

z <-seq(min(x),max(x),length.out=11)
lines(z,dea_est(x,b3,z,type="dea"),lty=1,col="red")
lines(z,dea_est(x,b3,z,type="fdh"),lty=2,col="blue")
lines(z,dea_est(x,b3,z,type="lfdh"),lty=3,col="green")
legend("topleft",legend=c("dea","fdh","lfdh"),col=c("red","blue","green"),lty=1:3)


• First of all thanks for your answer! Nonetheless, I am trying not to best fit/approximate the 1-ECDF overall but create only e.g. 5 blocks like i was trying to show in the picture. What I am saying is, that I am trying to find a best-possible fit of these 5 blocks to the function, not create a best-possible fit approximation overall. – BenR May 28 '14 at 14:19
• It's not clear from your question wether you also want to estimate the knot location. If not, that would be DEA with 5 knots. Again, the result could be pretty similar, except that DEA is implemented in C and thus presumably much faster to run. – user603 May 28 '14 at 14:22
• Thanks again! I actually would love to estimate the knot location if possible. – BenR May 28 '14 at 14:42
• estimating the knot location is certainly easier to do within the non parametric frontier regression framework (because it amount to adding a penalty to the objective function, exactly as is done in the usual, i.e. non-frontier non parametric regression framework) but to the best of my knowledge this is not widely implemented (yet) [we had the author of one of the packages I cite for a talk recently and he mentioned this, the adaptation to unknown knot location of the proposed method]. – user603 May 28 '14 at 16:46
• Again, this is very helpful. I think you are absolutely right. In questions sometimes less seems to be more... I will try to remodel my question although your answer seems to fit my problem. Somehow, I don't fully understand how to set the number of knots in the DEA. In your example it surely relates to the z?! – BenR May 29 '14 at 9:46