[Python] Is there (a Python implementation of) an average-consistent interpolation method?

I have a time series of values at low temporal resolution where the values represent an average of the respective surrounding intervals. Let's use hourly resolution for a simple example (in Python):

import pandas as pd
s_resample_test = pd.Series([0, 1, 2, 1, 0, -1, -2, -1], index=pd.date_range("2025-01-01 00:30:00", periods=8, freq="1h"))


So the value 0 at 00:30:00 represents an average for the interval 00:00:00 to 01:00:00, etc.

I wish to up-sample this sequence to a higher temporal resolution in such a way that the average of the up-sampled values is equal to the corresponding value from the original time series.
For example, when up-sampling the above time series, I want the average of the up-sampled values in the interval from 00:00:00 to 01:00:00 to be 0, the values in the interval from 01:00:00 to 02:00:00 to average to 1, etc.

Which interval do the points at the boundaries (00:00:00, 01:00:00, ...) belong to? I would consider splitting them between both intervals and half-weighting them in the average. To check solutions, I want my interpolated time series to obey (avg_mtx is the most precise way I can express my average constraint):

import numpy as np

init_mtx = np.eye(32)[0::4,:]
avg_mtx = np.empty((init_mtx.shape[0], init_mtx.shape[1] + 4))
for idx in range(init_mtx.shape[0]):
avg_mtx[idx, :] = np.convolve([.5, 1, 1, 1, .5], init_mtx[idx,:])
avg_mtx /= avg_mtx.sum(axis=1).reshape(-1,1)
avg_mtx = avg_mtx[1:-1,2:-5] # ignore boundary effects

s_resampled = desired_resampling_method(s_resample_test)

assert np.allclose(avg_mtx @ s_resampled, s_resample_test[1:-1])


For simplicity we can just say that the intervals are half-open: [00:00:00, 01:00:00), etc. This may facilitate a solution.

If I for example up-sample to quarterly resolution using:

s_resample_test.resample("15min").interpolate(method="polynomial", order=2)


or:

s_resample_test.resample("15min").interpolate(method="linear")


the result does not generally obey this average property. Neither does Fourier-based resampling using scipy.signal.resample.

Does an algorithm for the kind of interpolation I am looking for exist? And is there an implementation of it in Python?

A piecewise constant solution is probably the simplest to achieve, but I would like something smoother if possible. For example, I am speculating whether this can be formulated as a constrained smoothing spline approximation problem.

If not, I would also appreciate suggestions for a re-normalisation approach. Also, let us not necessarily worry about edge effects for now.

Continue reading...