ENH: Actually use rolling windows for EWMA in MACD

tests/pipeline/test_technical.py

@@ -5,6 +5,7 @@
 import numpy as np
 import pandas as pd
 import talib
+from numpy.random import random_integers
 
 from zipline.lib.adjusted_array import AdjustedArray
 from zipline.pipeline.data import USEquityPricing
@@ -16,6 +17,7 @@
     LinearWeightedMovingAverage,
     RateOfChangePercentage,
     TrueRange,
+    MovingAverageConvergenceDivergence
 )
 from zipline.testing import parameter_space
 from zipline.testing.fixtures import ZiplineTestCase
@@ -403,3 +405,78 @@ def test_tr_basic(self):
 
         tr.compute(today, assets, out, highs, lows, closes)
         assert_equal(out, np.full((3,), 2.))
+
+
+class MovingAverageConvergenceDivergenceCase(ZiplineTestCase):
+    def test_MACD_window_length_generation(self):
+        signal_period = random_integers(1, 90)
+        fast_period = random_integers(signal_period+1, signal_period+100)
+        slow_period = random_integers(fast_period+1, fast_period+100)
+        ewma = MovingAverageConvergenceDivergence(
+            fast_period=fast_period,
+            slow_period=slow_period,
+            signal_period=signal_period,
+        )
+        assert_equal(
+            ewma.window_length,
+            slow_period+signal_period-1,
+        )
+
+    def test_moving_average_convergence_divergence(self):
+        fast_period = 3
+        slow_period = 8
+        signal_period = 2
+
+        macd = MovingAverageConvergenceDivergence(
+            fast_period=fast_period,
+            slow_period=slow_period,
+            signal_period=signal_period,
+        )
+
+        today = pd.Timestamp('2016', tz='utc')
+        nassets = macd.window_length
+        assets = pd.Index(np.arange(nassets))
+        days_col = np.arange(start=-.05,
+                             stop=.01*nassets-.05,
+                             step=.01)[:, np.newaxis]
+        close = np.logspace(start=.01, stop=.10, num=nassets) - 1 + days_col
+
+        dtype = [
+            ('macd', 'f8'),
+            ('signal', 'f8'),
+            ('hist', 'f8'),
+        ]
+        out = np.recarray(
+            shape=(nassets,),
+            dtype=dtype,
+            buf=np.empty(shape=(nassets,), dtype=dtype),
+        )
+        macd.compute(
+            today,
+            assets,
+            out,
+            close,
+            fast_period,
+            slow_period,
+            signal_period,
+        )
+
+        expected_macd = np.array([0.01691553] * nassets)
+        expected_signal = np.array([0.01691553] * nassets)
+        expected_hist = np.array([0] * nassets)
+
+        np.testing.assert_almost_equal(
+            out.macd,
+            expected_macd,
+            decimal=8
+        )
+        np.testing.assert_almost_equal(
+            out.signal,
+            expected_signal,
+            decimal=8
+        )
+        np.testing.assert_almost_equal(
+            out.hist,
+            expected_hist,
+            decimal=8
+        )

zipline/pipeline/factors/__init__.py

@@ -32,6 +32,8 @@
     TrueRange,
     VWAP,
     WeightedAverageValue,
+    MovingAverageConvergenceDivergence,
+    MACD,
 )
 
 __all__ = [
@@ -62,4 +64,6 @@
     'TrueRange',
     'VWAP',
     'WeightedAverageValue',
+    'MovingAverageConvergenceDivergence',
+    'MACD'
 ]

zipline/pipeline/factors/technical.py

@@ -23,6 +23,7 @@
     sum as np_sum,
     nan
 )
+from numpy.lib.stride_tricks import as_strided
 from numexpr import evaluate
 
 from zipline.pipeline.data import USEquityPricing
@@ -40,8 +41,6 @@
 )
 from .factor import CustomFactor
 
-from talib import MACD
-
 
 class Returns(CustomFactor):
     """
@@ -448,11 +447,6 @@ def compute(self, today, assets, out, data, decay_rate):
         out[:] = sqrt(variance * bias_correction)
 
 
-# Convenience aliases.
-EWMA = ExponentialWeightedMovingAverage
-EWMSTD = ExponentialWeightedMovingStdDev
-
-
 class BollingerBands(CustomFactor):
     """
     Bollinger Bands technical indicator.
@@ -688,8 +682,14 @@ def compute(self, today, assets, out, highs, lows, closes):
         )
 
 
+def rolling_window(a, window):
+    # used to calculate rolling statistics using np.array
+    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
+    strides = a.strides + (a.strides[-1],)
+    return as_strided(a, shape=shape, strides=strides)
 
-class MovingAverageConvergenceDivergence(CustomFactor):
+
+class MovingAverageConvergenceDivergence(_ExponentialWeightedFactor):
     """
     Moving Average Convergence/Divergence (MACD)
     https://en.wikipedia.org/wiki/MACD
@@ -700,7 +700,8 @@ class MovingAverageConvergenceDivergence(CustomFactor):
     trend in a stock's price.
 
     **Default Inputs:** :data:`zipline.pipeline.data.USEquityPricing.close`
-    **Default Window Length:** None
+    **Default Window Length:** Window length is automatically calculated as the
+    sum of slow_period and signal_period.
 
     Parameters
     ----------
@@ -714,7 +715,7 @@ class MovingAverageConvergenceDivergence(CustomFactor):
     Returns
     -------
     MACD: The difference between "fast" EMA and "slow" EMA.
-    signal: The signal_period length period EMA of the MACD line.
+    signal: The EMA of the MACD line using `signal_period` as span.
     hist: Difference between MACD and signal. (Divergence series)
     """
     inputs = [USEquityPricing.close]
@@ -728,32 +729,53 @@ def __new__(cls,
                 *args,
                 **kwargs):
         return super(MovingAverageConvergenceDivergence, cls).__new__(
+            cls,
             fast_period=fast_period,
             slow_period=slow_period,
             signal_period=signal_period,
-            window_length=slow_period + signal_period,
+            window_length=slow_period + signal_period - 1,
             *args, **kwargs
         )
 
-    def calculate_macd(self, col, fast, slow, signal):
+    def calculate_ewa(self, data, length):
+        decay_rate = 1.0 - (2.0 / (1.0 + length))
+        return average(data,
+                       axis=1,
+                       weights=self.weights(length, decay_rate))
+
+    def rolling_window(self, a, window):
+        # used to calculate rolling statistics using np.array
+        shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
+        strides = a.strides + (a.strides[-1],)
+        return as_strided(a, shape=shape, strides=strides)
+
+    def calculate_macd(self, col):
         try:
-            macd, sig, hist = MACD(col,
-                                   fastperiod=fast,
-                                   slowperiod=slow,
-                                   signalperiod=signal)
-            return macd[-1], sig[-1], hist[-1]
+            slow_EWA = self.calculate_ewa(
+                self.rolling_window(
+                    col,
+                    self.params['slow_period']
+                ),
+                self.params['slow_period'])
+            fast_EWA = self.calculate_ewa(
+                self.rolling_window(
+                    col,
+                    self.params['fast_period']
+                )[-self.params['signal_period']:],
+                self.params['fast_period'])
+            macd = fast_EWA - slow_EWA
+            signal_line = self.calculate_ewa(
+                macd.reshape(-1, self.params['signal_period']),
+                self.params['signal_period'])
+            hist = macd[-1] - signal_line
+            return macd[-1], signal_line[-1], hist[-1]
         except:
             return nan, nan, nan
 
     def compute(self, today, assets, out, close, fast_period, slow_period,
                 signal_period):
-        n = len(close)
-        macd, sig, hist = zip(*map(self.calculate_macd,
-                                   close.T,
-                                   [fast_period]*n,
-                                   [slow_period]*n,
-                                   [signal_period]*n))
-        out.MACD[:] = macd
+        macd, sig, hist = zip(*map(self.calculate_macd, close.T))
+        out.macd[:] = macd
         out.signal[:] = sig
         out.hist[:] = hist
 
@@ -770,5 +792,12 @@ class AnnualVolatility(CustomFactor):
     """
     inputs = [USEquityPricing.close]
 
-    def compute(self, today, assets, out, closes):
-        out[:] = nanstd(closes, ddof=1, axis=0) * (252 ** 0.5)
+    params = {'annualization_factor': 252}
+
+    def compute(self, today, assets, out, closes, annualization_factor):
+        out[:] = nanstd(closes, ddof=1, axis=0) * (annualization_factor ** 0.5)
+
+# Convenience aliases.
+EWMA = ExponentialWeightedMovingAverage
+EWMSTD = ExponentialWeightedMovingStdDev
+MACD = MovingAverageConvergenceDivergence