BUG: Possible bug in rolling standard deviation calculation when using custom weights (pandas.core.window.rolling.Window.std)

[jevkus]

Pandas version checks

• I have checked that this issue has not already been reported.

• I have confirmed this bug exists on the latest version of pandas.

• I have confirmed this bug exists on the main branch of pandas.

Reproducible Example

import pandas as pd
import numpy as np

np.random.seed(42)
df = pd.Series(np.random.normal(size=1200))
df_std_stock = (
    df.rolling(window=220, min_periods=100, win_type='exponential')
      .std(tau=-(25 / np.log(2)), center=0, sym=False)
)

Issue Description

Hey, Team! The series of standard deviations produced by the rolling calculation with custom weights does not look right and I hope that a visual demonstration illustrates that best.

The blue line is the result of the calculation above and it is visible how after window observations the series no longer produces sensible results. Odd and strong changes in behavior are visible every pass over a window multiple, which is even more suspicious. The orange series is what I would expect the result to be and the calculation snippet that produced it is in the box below. The algorithm for the expected behavior is known to be problematic though.

Thank you for taking the time to consider my report and please let me know if further details are required. Lover your package and really hope that his helps making it even better.

Expected Behavior

import pandas as pd
import numpy as np

np.random.seed(42)
df = pd.Series(np.random.normal(size=1200))
df_std_custom = (
    (df ** 2).rolling(window=220, min_periods=100, win_type='exponential')
             .mean(tau=-(25 / np.log(2)), center=0, sym=False)
)
df_std_custom_mean = (
    df.rolling(window=220, min_periods=100, win_type='exponential')
      .mean(tau=-(25 / np.log(2)), center=0, sym=False)
)
df_std_custom = (df_std_custom - df_std_custom_mean**2).pow(0.5)

Installed Versions

INSTALLED VERSIONS

commit : 37ea63d
python : 3.9.16.final.0
python-bits : 64
OS : Windows
OS-release : 10
Version : 10.0.19044
machine : AMD64
processor : Intel64 Family 6 Model 141 Stepping 1, GenuineIntel
byteorder : little
LC_ALL : None
LANG : en_US.UTF-8
LOCALE : English_United States.1252

pandas : 2.0.1
numpy : 1.22.4
pytz : 2023.3
dateutil : 2.8.2
setuptools : 66.0.0
pip : 23.0.1
Cython : 0.29.34
pytest : 7.3.1
hypothesis : 6.75.3
sphinx : None
blosc : None
feather : None
xlsxwriter : None
lxml.etree : None
html5lib : None
pymysql : None
psycopg2 : None
jinja2 : 3.1.2
IPython : 8.13.2
pandas_datareader: None
bs4 : 4.12.2
bottleneck : None
brotli : None
fastparquet : None
fsspec : None
gcsfs : None
matplotlib : 3.7.1
numba : None
numexpr : 2.8.4
odfpy : None
openpyxl : None
pandas_gbq : None
pyarrow : 12.0.0
pyreadstat : None
pyxlsb : None
s3fs : None
scipy : 1.7.1
snappy : None
sqlalchemy : None
tables : 3.8.0
tabulate : None
xarray : None
xlrd : None
zstandard : None
tzdata : 2023.3
qtpy : None
pyqt5 : None

[topper-123]

This looks like a duplicate of #53289. Can you confirm?

[jevkus]

Thank you for the question @topper-123. I think these are separate issues. If I use my example with pd.rolling.std as in there, I get a match.

import pandas as pd
import numpy as np

np.random.seed(42)
df = pd.Series(np.random.normal(size=1200))

df_std = df.rolling(window=220).std(ddof=0)
df_std_np = df.rolling(window=220).apply(np.std)
df_custom = (df ** 2).rolling(window=220).mean()
df_custom = (df_custom - df.rolling(window=220).mean()**2).pow(0.5)

The linked issues suggests a problem with data type as it involves floats going over a wide range of magnitudes. Also, #53289 does not involve custom weights and possibly different code paths are involved. In the latter, code variance is calculated with this

pandas/pandas/_libs/window/aggregations.pyx

Lines 415 to 416 in 37ea63d

	def roll_var(const float64_t[:] values, ndarray[int64_t] start,
	ndarray[int64_t] end, int64_t minp, int ddof=1) -> np.ndarray:

while in this issue variance is calculated here

pandas/pandas/_libs/window/aggregations.pyx

Lines 1667 to 1668 in 37ea63d

	def roll_weighted_var(const float64_t[:] values, const float64_t[:] weights,
	int64_t minp, unsigned int ddof):

[topper-123]

Ok, thanks for looking into it.

[jevkus]

Hey @WillAyd, @mroeschke and @jreback! I hope not to bother and I am sorry for being so late to spot this, but why was PR #55153 closed? Is there a chance this could be reopen? I think @ihsansecer made a very compelling case for the fix, which may have been missed because of the fragmented nature of PR comments.

It is true that the new algorithm scales differently than the current one. However, the current algorithm is incorrect for rolling variance estimates and cannot produce meaningful results, so it is hard to say that this is a regression. @ihsansecer also points out that it is unlikely that there is a better solution, and I have to agree after reviewing standard references online. Most standard references, such as West (1979), are meant for estimate updates where the observed data/weight tuples do not change and new tuples are coming in, which is not the case in rolling calculations. In a rolling calculations with arbitrary weights, one needs to reweigh every observation once a window moves, so some nested looping needs to take place.

In addition, scaling of the new algorithm, which is O(w*N) where w is the window size and N is the size of the dataset, can still be practical. We use the workaround that I shared in the original post in practice. It uses the V[X] = E[X^2] - E[X]^2 relationship and pandas.DataFrame.rolling.mean method, and the calculation of the mean has the same scaling as the new variance calculations would have as @ihsansecer points out. Given that this scales linearly with the window for a fixed dataset, it serves us well with datasets up to billions of entries, but this workaround is generally problematic.

Besides the issue of a more complex implementation of V[X] = E[X^2] - E[X]^2, the approach can be unstable due to cancellation of large terms, see this Wiki page for reference. We took these risks because standard deviations are typically larger than the means in our applications and we had no reasonable alternatives. This is still a risk and we are always mindful of it once new applications are coming up. Weighted rolling variance estimates are very common in several domains, so I suspect we are not alone in this position. This suboptimal implementation has to be used by others too if meaningful results are needed.

Thank you in advance for taking the time to consider these points and to respond. Let me know if I missed something important that could be helpful for making progress with this issue. I really appreciate the work of the community and am happy to help to the best of my abilities.

[kaixiongg]

Hey @jevkus, I've also noticed that PR #55153 was closed, and the explanation by @ihsansecer makes sense to me. It seems there isn't an O(n) online rolling method for weighted mean/variance. However, if there is some pattern in the weights, such as a linear decay given the weights [1, 2, 3, 4...], a faster O(n) method for the weighted mean might be possible, though it may not be numerically stable.

[jevkus]

@kaixiongg, the explanation makes sense to me too and exponential weights in my example would fit this case. General applications of rolling, however, cannot assume that and additional complexity is perhaps OK to get to a correct result.