BUG: ewmvar/std/cov() de-biasing formula incorrect

[seth-p]

ewmvar() and ewmcov() contain the following code:

    if not bias:
        result *= (1.0 + 2.0 * com) / (2.0 * com)

Where does this scaling formula come from? Does the number of observations not enter in to the bias of the (sample) variance estimation? I would imagine that perhaps the formula above is for the case that have an infinite sample, but for a finite sample I would have thought that the number of observations would enter into things.

[seth-p]

Indeed, it looks like the existing scaling factor (1.0 + 2.0 * com) / (2.0 * com) is the special case of the "V1^2 / (V1^2 - V2)" formula of http://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Weighted_sample_variance -- or equivalently, the formula under gsl_stats_wvariance at http://www.gnu.org/software/gsl/manual/html_node/Weighted-Samples.html -- when have an infinite number of data points.

I think the code should be updated to calculate the correct formula for a finite sample. Though I think doing so may require writing a separate algos.ewmvar() analogous to the existing algos.roll_var() (as opposed to simply calculating ewmvar(x) = ewma(x*x) - ewma(x)**2).

[jreback]

@seth-p no idea why/when this came to be. If you can figure that out you are 90% there.

[seth-p]

Yep, numerical testing in Excel shows that the "V1^2 / (V1^2 - V2)" formula approaches (1.0 + 2.0 * com) / (2.0 * com) as the number of samples increases -- assuming adjust=True and there are no missing values.

I think I can implement the "V1^2 / (V1^2 - V2)" bias correction formula (which will depend on adjust, ignore_na, and the pattern of missing values) without resorting to writing a new algos.ewmvar(). So shouldn't be complicated.

[seth-p]

CC'ing @jaimefrio, in case of interest.