# Resampling/Interpolating monthly rates to daily rate estimates in R

I'm really not sure what to search for. If the answer to this is googleable, I'd be happy to hear what I should google.

I have a dataset of energy meter data. Readings are taken at roughly monthly intervals, though it is irregular. The data describe total energy use since the previous reading was taken. I want to resample that data to estimate usage for each month (the first to the last day). I'm looking into linear and spline interpolation, but both methods lead to high error rates. As shown below, integrating either of these methods for the known energy use of an interval results in some very large differences.

Is there a way to generate an interpolation for data like this while still being able to match the function over known ranges?

# Create Sample Meter Timeseries
ts <- data.frame(start_date = as.Date(c("2015-09-18",
"2015-10-20",
"2015-11-25",
"2015-12-23",
"2016-01-22")),
end_date = as.Date(c("2015-10-20",
"2015-11-25",
"2015-12-23",
"2016-01-22",
"2016-02-21")),
energy_use_kWh = c(3211,172,5566,8120,1344))

# Convert interval use into daily rate
ts$rate <- ts$energy_use_kWh/as.numeric(ts$end_date - ts$start_date)

# Generate midpoint of each interval
ts$mid_date <- ts$start_date + (ts$end_date - ts$start_date)/2

# Generate linear and spline interpolations between interval midpoints
approxfun_fun <- approxfun(x = ts$mid_date, y = ts$rate, rule = 2)
splinefun_fun <- splinefun(x = ts$mid_date, y = ts$rate, method = 'natural')

# Compare
ts$approx_estimate <- apply(ts[,c('start_date','end_date')], MARGIN = 1, function(x) integrate(f = approxfun_fun, lower = as.Date(x[1]), upper = as.Date(x[2]))$value)

ts$spline_estimate <- apply(ts[,c('start_date','end_date')], MARGIN = 1, function(x) integrate(f = splinefun_fun, lower = as.Date(x[1]), upper = as.Date(x[2]))$value)

result <- ts[c("start_date","end_date","energy_use_kWh","approx_estimate","spline_estimate")]
result$approx_error <- with(result, abs(energy_use_kWh - approx_estimate)/energy_use_kWh) result$spline_error <- with(result, abs(energy_use_kWh - spline_estimate)/energy_use_kWh)


EDIT: I have implemented GeoMatt22's suggestion over my sample data and found some interesting results...

Taking the integral over the known billing ranges matches perfectly with the known data, so that is definitely a success. The cumulative energy use graph looks pretty reasonable, though the derivative graph definitely looks a little strange. Within a billing period, the rate change is smooth and continuous, but there are some cases of instantaneous usage rate changes, where one day it is decreasing rapidly and the next day it is increasing rapidly.

This is the best answer so far, even though the day-to-day rates look a little funky. I'd welcome improvements that address this, but I am also happy with what I've got so far.

Edit 2: I tried changing the spline method from "monoH.FC", which computes a monotone Hermite spline, to "Hyman", which computes a monotone cubic spline using Hyman filtering. The results I got are a bit more continuous, and pass my personal eye-test a little better, though it's still not nicely continuous.

Edit 3: I built a function for the C2 monotone interpolant, as suggested by GeoMatt22. It took me a long time, but I got it to work! Method is from C2 rational quadratic spline interpolation to monotonic data (1983).

If you want more info on the code, I'd be happy to share. Currently in development, but it seems to be working well. The derivative curve, which corresponds to the daily use rate, has no unsightly kinks and integrating that function over each billing interval reproduces the known values exactly.

• In order for the end of month usage to agree with the daily usage, the average daily usage times the number of days must be set equal to the monthly usage. Thus, the simplest model would be to set the daily usage to the monthly usage divided by the number of days in that month. This would produce a step function, but, it would also conserve usage. One can smooth the step function using a stationary smooth like 1/4,1/2,1/4, or 1/16,1/4,3/8,1/4,1/16 (i.e., 1,2,1 divided by 4, or 1,4,6,4,1 divided by 16), but that will not conserve monthly totals exactly.
– Carl
Sep 26, 2016 at 19:35
• This was the first option I tried, since it seemed straightforward and conservative. I couldn't come up with a compact way of doing it though, so I didn't include it in my examples. The main downside I found with it is that usage probably isn't consistent over a billing interval. If there's a low November bill, a high December bill, and an even higher January bill, then it would probably be safe to assume that usage towards the end of the December bill is higher than usage at the beginning of the bill (i.e. the portion of the bill in November). Sep 26, 2016 at 19:45
• Personally, I would just write code that conserves energy (literally) and interpolates. For example, starting at day one and connecting to last month's energy usage, draw a triangle to the middle of the month, and from the apex of that triangle to the end of the month such that the area under the curve during that month is the energy usage for that month. If you cannot fathom the geometry, let me know, it is just simple algebra.
– Carl
Sep 26, 2016 at 19:55
• By day one do you mean the first of the month or the start date of the billing interval? The monthly energy usage is not known. Also, by triangle, do you mean that the first point is at zero? I suppose I cannot fathom the geometry. Sep 26, 2016 at 20:19
• My issues with integrate were probably due to my lack of any R skills whatsoever :) As I noted in my answer, the only constraint on your problem is that the $E(t)$ function must be non-decreasing, $E'(t)\geq 0$. This can be accomplished by any monotone interpolant. A standard cubic spline has continuous 1st and second derivatives, and is the smoothest possible piecewise cubic interpolant of the data. However this causes the common "overshoot" problems seen in e.g. the Wikipedia figure. The monotone version prevents this, but has a discontinuous 2nd derivative. This is why you see those kinks. Sep 27, 2016 at 19:42

If I understand the question correctly, the idea is to resample the energy-usage rates conservatively. To ensure conservative resampling, you should resample the extensive quantity ("mass" = cumulative energy used) rather than the intensive quantity ("density" = usage rate). This is very similar to how resampling a probability density is tricky, but resampling a cumulative distribution is straightforward (i.e. no coordinate-change adjustment is required).

In the current case, we have a time series of cumulative energy usage $(E_k,t_k)$. The original question is phrased in terms of the average energy-usage rate (power), which is the ratio of first differences, i.e. $\bar{r}_k=\Delta E_k/\Delta t_k$, and notes difficulty in conservative resampling of the $\bar{r}_k$ time series.

However, if we consider the energy series itself, then we have $$\Delta E_k = E_{k+1}-E_k = \int_{t_k}^{t_{k+1}}r[t]dt = \bar{r}_k\Delta t_k$$

The energy series itself $E(t)$ can be interpolated using any monotone scheme (e.g. linear interpolation, or monotone cubic). The monotone requirement ensures that it will always be non-decreasing through time: $E_k\leq E_{k+\phi}\leq E_{k+1}$ for $\phi\in[0,1]$.

Once this is done, the new energy series can be differenced to get average usage rates over the new time intervals.

Summary: the average usage rate is by definition $\bar{r}=\Delta E/\Delta t$, so if you interpolate the cumulative energy $E(t)$ (monotonically) then you will automatically have conservative results.

I do not use R, but skimming the help, it looks like you can do something like:

nt <- length(ts$start_date) t <- c(ts$start_date,ts$end_date[nt]) E <- cumsum(c(0,ts$energy_use))
Espline <- splinefun(x = t, y = E, method = 'monoH.FC')
dEdt_spline <- function(t) Espline( t , deriv = 1 )


Then you can evaluate the average power consumption as $\langle r\rangle_{t\in[t_1,t_2]} = \frac{E(t_2)-E(t_1)}{t_2-t_1}$, and the "instantaneous" power consumption with $E'(t)$.

(Note: I quickly tried this on R-fiddle and it seemed to work, but your integrate test still did not work. I strongly believe this must be due to some code error, either on my part or in the R libraries. That is, by design $E(t)$ has been fit with a monotone interpolating spline, which has an analytic derivative, that itself has an analytic integral, as they are piecewise polynomial functions. Most likely the inconsistency is due to my lack of R knowledge, or it could be due to numerical approximations used in the spline calculus functions.)

Update: As expected, the above was due to my lack of R knowledge, as shown in updated question. (I had literally never written any R before Googling to do the above, so not too shabby!) Note also that as seen there, the monotone cubic spline functions will have a discontinuous second derivative (seen as kinks in the $E'(t)$ plots). This could be avoided with a monotone C2 interpolant (e.g. this), though I do not know what R package this might be in.

Note on monotone interpolation: This simply means that the interpolation does not introduce any new local maxima/minima not already present in the data. For example the following picture from Wikipedia demonstrates how the standard cubic spline is not monotone

Public Domain, https://en.wikipedia.org/w/index.php?curid=9051137

• The existence of this SO question might suggest that R's integrate() is not smart enough to treat piecewise polynomials analytically. (Not sure if there is an alternative?) Regardless, the average rates are easily obtained from the spline of the cumulative energy directly, so integration is not needed (and is more roundabout anyway!). Sep 27, 2016 at 16:25
• I actually had good success with your function and integrate() - the integrations matched perfectly with the known ranges. Sep 27, 2016 at 19:10
• @jmartindill great! If you found my answer useful, you can upvote and/or accept it. Sep 27, 2016 at 19:43
• Alright I accepted it. Check out the edits I made to my main post, if you're interested. I implemented your method, but tried another spline method that uses a cubic spline with Hyman filtering. Sep 27, 2016 at 19:50
• You would need a "monotone C2 interpolant" to get rid of the kinks. I do not know what solutions may exist in R. A quick Google suggests the method in C2 rational quadratic spline interpolation to monotonic data (1983) would be appropriate (e.g. as discussed here). You may have to implement it yourself. (Or ask about "monotone C2 interpolant" on an R list?) Sep 27, 2016 at 20:00

What I did is calculate the usage at mid month needed to balance the end month kw usage calculation using linear interpolation. This works out to be a simple formula: $kw_{m-1/2}=\frac{3}{2}kw_m-\frac{1}{2}kw_{m-1}$. That is, the mid month consumption is 3/2 of the end month consumption minus 1/2 of the prior month's utilization. Why? Simple geometry. Note that the area of each trapezoid under the curve is the average height times the base. Note that the area of the two half- month trapezoids culminating in the end of month usage is thus

$\frac{1}{2}\frac{kw_{m-1}+\frac{3}{2}kw_m-\frac{1}{2}kw_{m-1}}{2}+\frac{1}{2}\frac{\frac{3}{2}kw_m-\frac{1}{2}kw_{m-1}+kw_m}{2}=kw_m.$

To show this, I made up pseudo data because I didn't have any. Here are the step by step calculations.

Pseudo-data and calculations:

month end   kw  days    kw/day       mid month          mid month usage
01/10/2017  546                                         3/2*kw(m)-1/2*kw(m-1)
01/11/2017  294 31     17.61290323   16/10/2017 12:00   168
01/12/2017  493 30     9.8           16/11/2017 0:00    592.5
01/01/2018  593 31     15.90322581   16/12/2017 12:00   643


Putting together the true usages and mid month usages in a single table give us a table of all values:

dates               kw
01/10/2017          546
16/10/2017 12:00    168
01/11/2017          294
16/11/2017 0:00     592.5
01/12/2017          493
16/12/2017 12:00    643
01/01/2018          593


We then proceed to plot that table:

Note, we could further require that the slope for each month and mid month be the same approached from the left or right, or we could solve the same problem a bazillion other ways. No matter how it is solved, there will be "slumps" or "bumps" in between the months to account, respectively, for the greater or lesser utilization during the prior month. For real life utilization, the actual curve is closer to a step function with many mostly tiny steps and some big ones too, so there is no way we can emulate it exactly.

• Carl, this seems like another ill-posed deconvolution problem. I have to look more closely to be sure, but it looks like a problem I have done before: If you have a smoother $A$ such that $(Ax)_i=(x_i+x_{i+1})/2$, then it has a null-space of $x_k=(-1)^k$ (i.e. an oscillatory "zero-energy mode" that gives 0 average, though it is non-zero itself). This would explain the non-monotone interpolant you get. (This could be regularized with a Tikhonov matrix of $(Dx)_i=x_{i+1}-x_i$.) Sep 26, 2016 at 22:29
• @GeoMatt22 It occurred to me that there are many many ways of approaching this problem. However, what the OP needed was something simple as a take home message, so, I went with the second simplest possible solution that fulfills the stated requirements, where the simplest; a step function, would not be continuous thus not interpolated at all.
– Carl
Sep 26, 2016 at 23:30
• actually looking at the problem, I think the simplest way to go is to interpolate the cumulative energy usage time series, which will ensure conservative rates (so long as monotone interpolation is used, i.e. $E_t\leq E_{t+\phi} \leq E_{t+1}$ for $\phi \in (0,1)$.) Sep 27, 2016 at 2:09
• @GeoMatt22 Yup, but that not what the OP asked for. He wanted interpolation for data like this while still being able to match the function over known ranges. I only solved that for month to month. Now the question you have is "Why do that?" I think it is for accounting purposes. Now vote for me somewhere. I'm stuck at 980 and I want to see the +/- voting :p
– Carl
Sep 27, 2016 at 2:27
• I do not see how "match the function over known ranges" differs from my answer? From the OP's description, and their demo, the data they have is measurements of $E(t)$ to begin with. The rates are then derived from these. I will check around your posts :) Sep 27, 2016 at 2:35