I think it's time I closed this Q&A. After a couple of years since asking this question a possible solution came to my mind, which I elaborated into my master thesis. I actually wanted to publish it in some scientific journal, and link the article here, but that turned out to be more difficult than I thought and also not worth the effort since I decided to not pursue a career in academy. Anyhow, my work is not available anywhere, but the end product is here: https://github.com/c-foschi/mapinterval.
Some notes about my solution (from now on Mapinterval) and its relation to other answers:
- Mapinterval corresponds to @ElonPlotkin's answer when the number of computed nodes tend to infinity. In other words, Mapinterval finds the the function that minimizes the integral of the square first derivative, among all functions which first derivative exists on the whole domain.
- I find some similarities between my solution and @RobertDodier's. Mapinterval is also a spline - a quadratic one though, and is also a natural spline since its first derivative is 0 at both ends of the time span.
- Mapinterval uses spline computation methods, and it's fast. It's complexity is O(n).
- I once found a presentation somewhere where they found a solution for this same problem that I think corresponds completely to Mapinterval: they computed the cumulated time series, then they fit a natural cubic spline to that series, and then they computed the derivative of that function.
- One property that was central to my thesis was that Mapinterval is the maximum a-posteriori estimator for the series generating process if we assume that process was Wiener (Brownian Motion). In my thesis I actually wrote that it was the maximum likelihood estimator, but my professor said that my model was properly Bayesian instead. Anyway, this allow us to draw confidence/credibility bands around the estimate, which looks really cool:
