ENH: Add support for calculating EWMA with a time component (#34839)

Author: mroeschke

asv_bench/benchmarks/rolling.py

@@ -91,11 +91,18 @@ class EWMMethods:
     def setup(self, constructor, window, dtype, method):
         N = 10 ** 5
         arr = (100 * np.random.random(N)).astype(dtype)
+        times = pd.date_range("1900", periods=N, freq="23s")
         self.ewm = getattr(pd, constructor)(arr).ewm(halflife=window)
+        self.ewm_times = getattr(pd, constructor)(arr).ewm(
+            halflife="1 Day", times=times
+        )
 
     def time_ewm(self, constructor, window, dtype, method):
         getattr(self.ewm, method)()
 
+    def time_ewm_times(self, constructor, window, dtype, method):
+        self.ewm.mean()
+
 
 class VariableWindowMethods(Methods):
     params = (

doc/source/user_guide/computation.rst

@@ -1095,6 +1095,25 @@ and **alpha** to the EW functions:
   one half.
 * **Alpha** specifies the smoothing factor directly.
 
+.. versionadded:: 1.1.0
+
+You can also specify ``halflife`` in terms of a timedelta convertible unit to specify the amount of
+time it takes for an observation to decay to half its value when also specifying a sequence
+of ``times``.
+
+.. ipython:: python
+
+    df = pd.DataFrame({'B': [0, 1, 2, np.nan, 4]})
+    df
+    times = ['2020-01-01', '2020-01-03', '2020-01-10', '2020-01-15', '2020-01-17']
+    df.ewm(halflife='4 days', times=pd.DatetimeIndex(times)).mean()
+
+The following formula is used to compute exponentially weighted mean with an input vector of times:
+
+.. math::
+
+    y_t = \frac{\sum_{i=0}^t 0.5^\frac{t_{t} - t_{i}}{\lambda} x_{t-i}}{0.5^\frac{t_{t} - t_{i}}{\lambda}},
+
 Here is an example for a univariate time series:
 
 .. ipython:: python

doc/source/whatsnew/v1.1.0.rst

@@ -329,6 +329,7 @@ Other enhancements
 - :meth:`DataFrame.to_excel` can now also write OpenOffice spreadsheet (.ods) files (:issue:`27222`)
 - :meth:`~Series.explode` now accepts ``ignore_index`` to reset the index, similarly to :meth:`pd.concat` or :meth:`DataFrame.sort_values` (:issue:`34932`).
 - :meth:`read_csv` now accepts string values like "0", "0.0", "1", "1.0" as convertible to the nullable boolean dtype (:issue:`34859`)
+- :class:`pandas.core.window.ExponentialMovingWindow` now supports a ``times`` argument that allows ``mean`` to be calculated with observations spaced by the timestamps in ``times`` (:issue:`34839`)
 
 .. ---------------------------------------------------------------------------
 

pandas/_libs/window/aggregations.pyx

@@ -8,7 +8,7 @@ from libc.stdlib cimport malloc, free
 
 import numpy as np
 cimport numpy as cnp
-from numpy cimport ndarray, int64_t, float64_t, float32_t
+from numpy cimport ndarray, int64_t, float64_t, float32_t, uint8_t
 cnp.import_array()
 
 
@@ -1752,6 +1752,51 @@ def roll_weighted_var(float64_t[:] values, float64_t[:] weights,
 # ----------------------------------------------------------------------
 # Exponentially weighted moving average
 
+def ewma_time(ndarray[float64_t] vals, int minp, ndarray[int64_t] times,
+              int64_t halflife):
+    """
+    Compute exponentially-weighted moving average using halflife and time
+    distances.
+
+    Parameters
+    ----------
+    vals : ndarray[float_64]
+    minp : int
+    times : ndarray[int64]
+    halflife : int64
+
+    Returns
+    -------
+    ndarray
+    """
+    cdef:
+        Py_ssize_t i, num_not_nan = 0, N = len(vals)
+        bint is_not_nan
+        float64_t last_result
+        ndarray[uint8_t] mask = np.zeros(N, dtype=np.uint8)
+        ndarray[float64_t] weights, observations, output = np.empty(N, dtype=np.float64)
+
+    if N == 0:
+        return output
+
+    last_result = vals[0]
+
+    for i in range(N):
+        is_not_nan = vals[i] == vals[i]
+        num_not_nan += is_not_nan
+        if is_not_nan:
+            mask[i] = 1
+            weights = 0.5 ** ((times[i] - times[mask.view(np.bool_)]) / halflife)
+            observations = vals[mask.view(np.bool_)]
+            last_result = np.sum(weights * observations) / np.sum(weights)
+
+        if num_not_nan >= minp:
+            output[i] = last_result
+        else:
+            output[i] = NaN
+
+    return output
+
 
 def ewma(float64_t[:] vals, float64_t com, bint adjust, bint ignore_na, int minp):
     """
@@ -1761,9 +1806,9 @@ def ewma(float64_t[:] vals, float64_t com, bint adjust, bint ignore_na, int minp
     ----------
     vals : ndarray (float64 type)
     com : float64
-    adjust: int
-    ignore_na: bool
-    minp: int
+    adjust : int
+    ignore_na : bool
+    minp : int
 
     Returns
     -------
@@ -1831,10 +1876,10 @@ def ewmcov(float64_t[:] input_x, float64_t[:] input_y,
     input_x : ndarray (float64 type)
     input_y : ndarray (float64 type)
     com : float64
-    adjust: int
-    ignore_na: bool
-    minp: int
-    bias: int
+    adjust : int
+    ignore_na : bool
+    minp : int
+    bias : int
 
     Returns
     -------

pandas/core/generic.py

@@ -10518,6 +10518,7 @@ def ewm(
             adjust=True,
             ignore_na=False,
             axis=0,
+            times=None,
         ):
             axis = self._get_axis_number(axis)
             return ExponentialMovingWindow(
@@ -10530,6 +10531,7 @@ def ewm(
                 adjust=adjust,
                 ignore_na=ignore_na,
                 axis=axis,
+                times=times,
             )
 
         cls.ewm = ewm

pandas/core/window/ewm.py

@@ -1,18 +1,21 @@
+import datetime
 from functools import partial
 from textwrap import dedent
 from typing import Optional, Union
 
 import numpy as np
 
+from pandas._libs.tslibs import Timedelta
 import pandas._libs.window.aggregations as window_aggregations
-from pandas._typing import FrameOrSeries
+from pandas._typing import FrameOrSeries, TimedeltaConvertibleTypes
 from pandas.compat.numpy import function as nv
 from pandas.util._decorators import Appender, Substitution, doc
 
+from pandas.core.dtypes.common import is_datetime64_ns_dtype
 from pandas.core.dtypes.generic import ABCDataFrame
 
 from pandas.core.base import DataError
-import pandas.core.common as com
+import pandas.core.common as common
 from pandas.core.window.common import _doc_template, _shared_docs, zsqrt
 from pandas.core.window.rolling import _flex_binary_moment, _Rolling
 
@@ -32,7 +35,7 @@ def get_center_of_mass(
     halflife: Optional[float],
     alpha: Optional[float],
 ) -> float:
-    valid_count = com.count_not_none(comass, span, halflife, alpha)
+    valid_count = common.count_not_none(comass, span, halflife, alpha)
     if valid_count > 1:
         raise ValueError("comass, span, halflife, and alpha are mutually exclusive")
 
@@ -76,10 +79,17 @@ class ExponentialMovingWindow(_Rolling):
     span : float, optional
         Specify decay in terms of span,
         :math:`\alpha = 2 / (span + 1)`, for :math:`span \geq 1`.
-    halflife : float, optional
+    halflife : float, str, timedelta, optional
         Specify decay in terms of half-life,
         :math:`\alpha = 1 - \exp\left(-\ln(2) / halflife\right)`, for
         :math:`halflife > 0`.
+
+        If ``times`` is specified, the time unit (str or timedelta) over which an
+        observation decays to half its value. Only applicable to ``mean()``
+        and halflife value will not apply to the other functions.
+
+        .. versionadded:: 1.1.0
+
     alpha : float, optional
         Specify smoothing factor :math:`\alpha` directly,
         :math:`0 < \alpha \leq 1`.
@@ -124,6 +134,18 @@ class ExponentialMovingWindow(_Rolling):
     axis : {0, 1}, default 0
         The axis to use. The value 0 identifies the rows, and 1
         identifies the columns.
+    times : str, np.ndarray, Series, default None
+
+        .. versionadded:: 1.1.0
+
+        Times corresponding to the observations. Must be monotonically increasing and
+        ``datetime64[ns]`` dtype.
+
+        If str, the name of the column in the DataFrame representing the times.
+
+        If 1-D array like, a sequence with the same shape as the observations.
+
+        Only applicable to ``mean()``.
 
     Returns
     -------
@@ -159,6 +181,17 @@ class ExponentialMovingWindow(_Rolling):
     2  1.615385
     3  1.615385
     4  3.670213
+
+    Specifying ``times`` with a timedelta ``halflife`` when computing mean.
+
+    >>> times = ['2020-01-01', '2020-01-03', '2020-01-10', '2020-01-15', '2020-01-17']
+    >>> df.ewm(halflife='4 days', times=pd.DatetimeIndex(times)).mean()
+              B
+    0  0.000000
+    1  0.585786
+    2  1.523889
+    3  1.523889
+    4  3.233686
     """
 
     _attributes = ["com", "min_periods", "adjust", "ignore_na", "axis"]
@@ -168,20 +201,49 @@ def __init__(
         obj,
         com: Optional[float] = None,
         span: Optional[float] = None,
-        halflife: Optional[float] = None,
+        halflife: Optional[Union[float, TimedeltaConvertibleTypes]] = None,
         alpha: Optional[float] = None,
         min_periods: int = 0,
         adjust: bool = True,
         ignore_na: bool = False,
         axis: int = 0,
+        times: Optional[Union[str, np.ndarray, FrameOrSeries]] = None,
     ):
+        self.com: Optional[float]
         self.obj = obj
-        self.com = get_center_of_mass(com, span, halflife, alpha)
         self.min_periods = max(int(min_periods), 1)
         self.adjust = adjust
         self.ignore_na = ignore_na
         self.axis = axis
         self.on = None
+        if times is not None:
+            if isinstance(times, str):
+                times = self._selected_obj[times]
+            if not is_datetime64_ns_dtype(times):
+                raise ValueError("times must be datetime64[ns] dtype.")
+            if len(times) != len(obj):
+                raise ValueError("times must be the same length as the object.")
+            if not isinstance(halflife, (str, datetime.timedelta)):
+                raise ValueError(
+                    "halflife must be a string or datetime.timedelta object"
+                )
+            self.times = np.asarray(times.astype(np.int64))
+            self.halflife = Timedelta(halflife).value
+            # Halflife is no longer applicable when calculating COM
+            # But allow COM to still be calculated if the user passes other decay args
+            if common.count_not_none(com, span, alpha) > 0:
+                self.com = get_center_of_mass(com, span, None, alpha)
+            else:
+                self.com = None
+        else:
+            if halflife is not None and isinstance(halflife, (str, datetime.timedelta)):
+                raise ValueError(
+                    "halflife can only be a timedelta convertible argument if "
+                    "times is not None."
+                )
+            self.times = None
+            self.halflife = None
+            self.com = get_center_of_mass(com, span, halflife, alpha)
 
     @property
     def _constructor(self):
@@ -277,14 +339,23 @@ def mean(self, *args, **kwargs):
             Arguments and keyword arguments to be passed into func.
         """
         nv.validate_window_func("mean", args, kwargs)
-        window_func = self._get_roll_func("ewma")
-        window_func = partial(
-            window_func,
-            com=self.com,
-            adjust=self.adjust,
-            ignore_na=self.ignore_na,
-            minp=self.min_periods,
-        )
+        if self.times is not None:
+            window_func = self._get_roll_func("ewma_time")
+            window_func = partial(
+                window_func,
+                minp=self.min_periods,
+                times=self.times,
+                halflife=self.halflife,
+            )
+        else:
+            window_func = self._get_roll_func("ewma")
+            window_func = partial(
+                window_func,
+                com=self.com,
+                adjust=self.adjust,
+                ignore_na=self.ignore_na,
+                minp=self.min_periods,
+            )
         return self._apply(window_func)
 
     @Substitution(name="ewm", func_name="std")

pandas/tests/window/conftest.py

@@ -1,4 +1,4 @@
-from datetime import datetime
+from datetime import datetime, timedelta
 
 import numpy as np
 from numpy.random import randn
@@ -302,3 +302,9 @@ def series():
 def which(request):
     """Turn parametrized which as fixture for series and frame"""
     return request.param
+
+
+@pytest.fixture(params=["1 day", timedelta(days=1)])
+def halflife_with_times(request):
+    """Halflife argument for EWM when times is specified."""
+    return request.param