How Random is the Market? Testing the Random Walk Hypothesis

A mainstay of academic research into the market is the Random Walk Hypothesis (RWH). This is the idea that market moves are random and follow a normal distribution that can be easily described using a concept borrowed from physics called Brownian Motion.

This makes the market mathematics manageable, but is it true? Is the market really random?

If it is, then there’s little point to trying to beat it. But if it isn’t, then there are repeatable patterns that can be algorithmically exploited.

Thankfully, the issue of randomness is very important for fields like cryptography, so it is well studied and there are statistical tests that we can apply to market data to investigate this.

We’re going to borrow a few standard tests for randomness and apply it to historical data to see just how randome the markets really are.

Measuring Market Randomness

There are a host of randomness tests that have been developed over the years which look at binary sequences to determine whether or not a random process was used to generate these values. Notably, we have test suites such as the Diehard TestsTestU01NIST tests and others that have been published over the years.

We could run a large battery of tests (maybe we’ll get to that in a future post)to test our market data, but for now, we’ll just select three tests to see how the RWH holds up: runs test, discrete Fourier Transform test, and the Binary Matrix Rank test from the NIST suite.

Runs Test

If the market truly is random, then we shouldn’t see any dependence on previous prices; the market being up today should have no impact on what it will do tomorrow (and vice versa).

The runs test can help us look this aspect of randomness. It works by looking at the total number of positive and negative streaks in a sequence and checking the lengths.

We’ll take our prices and make all positive price changes into 1s and negative changes into 0s, and keep this binary vector as X. We’ll set n as the number of observations we have (e.g. n = len(X)). Then, to implement the runs test, we take the following steps (adapted from section 2.3 of the NIST Statistical Test Suite):

1. Compute the proportion of 1s in the binary sequence:

\pi = \frac{\sum_j X_j}{n}

2. Check the value \pi against the frequency test. It passes if: \mid \pi - 1/2 \mid < \tau, where \tau = \frac{2}{\sqrt{n}}. If the frequency test is failed, then we can stop and we don’t have a random sequence and we’ll set our P-value to 0. If we pass, then we can continue to step 3.

3. Compute our test statistic V_n where:

V_n = \sum_{k=1}^{n-1} r(k) + 1

where r(k) = 0 if X_k = X_{k+1}, otherwise r(k) = 1. So if we have the sequence [0, 1, 0, 0, 0, 1, 1], then this becomes: V_n = (1 + 1 + 0 + 0 + 1 + 0) + 1 = 4

4. Compute our P-value where:

p = erfc\bigg( \frac{ \mid V_n - 2n \pi (1 - \pi) \mid}{2 \pi (1-\pi) \sqrt{2n}} \bigg)

Note that erfc is the complementary error function (given below). Thankfully, this is available in Python with scipy.special.erfc(z):

erfc(z)=\frac{2}{\sqrt{\pi}}\int_z^{\infty}e^{-u^2}du

With all of that, we can now use our P-value to determine whether or not our sequence is random. If our P-value is below our threshold (e.g. 5%), then we reject the null hypothesis, which means we have a non-random sequence on our hands.

import numpy as np
from scipy.special import erfc

def RunsTest(x):
  # Convert input to binary values
  X = np.where(x > 0, 1, 0)
  n = len(X)
  pi = X.sum() / n
  # Check frequency test
  tau = 2 / np.sqrt(n)
  if np.abs(pi - 0.5) >= tau:
    # Failed frequency test
    return 0
  r_k = X[1:] != X[:-1]
  V_n = r_k.sum() + 1
  num = np.abs(V_n - 2 * n * pi * (1 - pi))
  den = 2 * pi * (1 - pi) * np.sqrt(2 * n)
  return erfc(num / den)

The NIST documentation gives us some test data to check that our function is working properly, so let’s drop that into our function and see what happens.

# eps from NIST doc
eps = '110010010000111111011010101000100010000101101' + \ 
  '0001100001000110100110001001100011001100010100010111000'
x = np.array([int(i) for i in eps])

p = RunsTest(x)
H0 = p > 0.01
# NIST P-value = 0.500798
print("Runs Test\n"+"-"*78)
if H0:
  print(f"Fail to reject the Null Hypothesis (p={p:.3f}) -> random sequence")
else:
  print(f"Reject the Null Hypothesis (p={p:.3f}) -> non-random sequence.")
Runs Test
------------------------------------------------------------------------------
Fail to reject the Null Hypothesis (p=0.501) -> random sequence

We get the same P-value, so we can be confident that our implementation is correct. Note also that NIST recommends we have at least 100 samples in our data for this test to be valid (i.e. $n \geq 100$).

Discrete Fourier Transformation Test

Our next test is the Discrete Fourier Transformation (DFT) test.

This test computes a Fourier Transform on the data and looks at the peak heights. If their are too many high peaks, then it indicates we aren’t dealing with a random process. It would take us too far afield to dive into the specifics of Fourier Transforms, but check out this post if you’re interested to go deeper.

Let’s get to the NIST steps. We have data (x) and we need to set a threshold, which is usually 95% as inputs.

1. We need to convert our time-series x into a sequence of 1s and -1s for positive and negative deviations. This new sequence is called \hat{x}.

2. Apply discrete Fourier Transform (DFT) to \hat{x}:

\Rightarrow S = DFT(\hat{x})

3. Calculate M = modulus(S') = \left| S \right|, where S' is the first n/2 elements in S and the modulus yields the height of the peaks.

4. Compute the 95% peak height threshold value. If we are assuming randomness, then 95% of the values obtained from the test should not exceed T.

T = \sqrt{n\textrm{log}\frac{1}{0.05}}

5. Compute N_0 = \frac{0.95n}{2}, where N_0 is the theoretical number of peaks (95%) that are less than T (e.g. if n=10, then N_0 = \frac{10 \times 0.95}{2} = 4.75).

6. Compute the P-value using the erfc function:

P = erfc \bigg( \frac{\left| d \right|}{\sqrt{2}} \bigg)

Just like we did above, we’re going to compare our P-value to our reference level and see if we can reject the null hypothesis – that we have a random sequence – or not. Note too that it is recommended that we use at least 1,000 inputs (n \geq 1000) for this test.

def DFTTest(x, threshold=0.95):
  n = len(x)
  # Convert to binary values
  X = np.where(x > 0, 1, -1)
  # Apply DFT
  S = np.fft.fft(X)
  # Calculate Modulus
  M = np.abs(S[:int(n/2)])
  T = np.sqrt(n * np.log(1 / (1 - threshold)))
  N0 = threshold * n / 2
  N1 = len(np.where(M < T)[0])
  d = (N1 - N0) / np.sqrt(n * (1-threshold) * threshold / 4)
  # Compute P-value
  return erfc(np.abs(d) / np.sqrt(2))

NIST gives us some sample data to test our implementation here too.

# Test sequence from NIST
eps = '110010010000111111011010101000100010000101101000110000' + \
  '1000110100110001001100011001100010100010111000'
x = np.array([int(i) for i in eps])
p = DFTTest(x)
H0 = p > 0.01
print("DFT Test\n"+"-"*78)
if H0:
  print(f"Fail to reject the Null Hypothesis (p={p:.3f}) -> random sequence")
else:
  print(f"Reject the Null Hypothesis (p={p:.3f}) -> non-random sequence.")
DFT Test
------------------------------------------------------------------------------
Fail to reject the Null Hypothesis (p=0.646) -> random sequence

Same as the NIST documentation, we reject the null hypothesis.

Binary Matrix Rank Test

We’ll choose one last test out of the test suite – the Binary Matrix Rank Test.

Steps:

1. Divide the sequence into 32 by 32 blocks. We’ll have N total blocks to work with and discard any data that doesn’t fit nicely into our 32×32 blocks. Each block will be a matrix consisting of our ordered data. A quick example will help illustrate, say we have a set of 10, binary data points: X = [0, 0, 0, 1, 1, 0, 1, 0, 1, 0] and we have 2×2 matrices (to make it easy) instead of 32×32. We’ll divide this data into two blocks and discard two data points. So we have two blocks (B_1 and B_2) that now look like:

B_1 = \begin{bmatrix}0 & 0 \\0 & 1\end{bmatrix} B_2 = \begin{bmatrix}1 & 0 \\1 & 0\end{bmatrix}

2. We determine the rank of each binary matrix. If you’re not familiar with the procedure, check out this notebook here for a great explanation. In Python, we can simply use the np.linalg.matrix_rank() function to compute it quickly.

3. Now that we have the ranks, we’re going to count the number of full rank matrices (if we have 32×32 matrices, then a full rank matrix has a rank of 32) and call this number F_m. Then we’ll get the number of matrices with rank one less than full rank which will be F_{m-1}. We’ll use N to denote the total number of matrices we have.

4. Now, we compute the Chi-squared value for our data with the following equation:

\chi^2 = \frac{(F_m-0.2888N)^2}{0.2888N} + \frac{(F_{m-1} - 0.5776N)^2}{0.5776N} + \frac{(N - F_m - F_{m-1} - 0.1336N)^2}{0.1336N}
  1. Calculate the P-value using the Incomplete Gamma Function, Q\big(1, \frac{\chi^2}{2} \big):
P = Q \bigg(1, \frac{\chi^2}{2} \bigg) = \frac{1}{\Gamma(a)} \int_x^{\infty} t^{a-1} e^{-1} = e^{\frac{\chi^2}{2}}

Scipy makes this last bit easy with a simple function call to scipy.special.gammaincc().

Don’t be intimidated by this! It’s actually straightforward to implement.

from scipy.special import gammaincc

def binMatrixRankTest(x, M=32):
  X = np.where(x > 0, 1, 0)
  n = len(X)
  N = np.floor(n / M**2).astype(int)
  # Create blocks
  B = X[:N * M**2].reshape(N, M, M)
  ranks = np.array([np.linalg.matrix_rank(b) for b in B])
  F_m = len(np.where(ranks==M)[0])
  F_m1 = len(np.where(ranks==M - 1)[0])
  chi_sq = (F_m - 0.2888 * N) ** 2 / (0.2888 * N) \
    + (F_m1 - 0.5776 * N) ** 2 / (0.5776 * N) \
    + (N - F_m - F_m1 - 0.1336 * N) ** 2 / (0.1336 * N)
  return gammaincc(1, chi_sq / 2)

If our P-value is less than our threshold, then we have a non-random sequence. Let’s test it with the simple example given in the NIST documentation to ensure we implemented things correctly:

eps = '01011001001010101101'
X = np.array([int(i) for i in eps])
p = binMatrixRankTest(X, M=3)
H0 = p > 0.01
print("Binary Matrix Rank Test\n"+"-"*78)
if H0:
  print(f"Fail to reject the Null Hypothesis (p={p:.3f}) -> random sequence")
else:
  print(f"Reject the Null Hypothesis (p={p:.3f}) -> non-random sequence.")
Binary Matrix Rank Test
------------------------------------------------------------------------------
Fail to reject the Null Hypothesis (p=0.742) -> random sequence

And it works! Note that in this example, we have a much smaller data set, so we set M=3 for 9-element matrices. This test is also very data hungry. They recommend at least 38 matrices to test. If we’re using 32×32 matrices, then that means we’ll need 38x32x32 = 38,912 data points. That’s roughly 156 years of daily price data!

Only the oldest companies and commodities are going to have that kind of data available (and not likely for free). We’ll press on with this test anyway, but take the results with a grain of salt because we’re violating the data recommendations.

Testing the RWH

With our tests in place, we can get some actual market data and see how well the RWH holds up. To do this properly, we’re going to need a lot of data, so I picked out some indices with a long history, a few old and important commodities, some of the oldest stocks out there, a few currency pairs, and Bitcoin just because.

Data from:

  • Dow Jones
  • S&P 500
  • Gold
  • Oil
  • USD/GBP
  • BTC/USD

One thing to note as well, we want to also run this against a baseline. For each of these I’ll be benchmarking the results against NumPy’s binomial sampling algorithm, which should have a high-degree of randomness.

I relied only on free sources so you can replicate this too, but more and better data is going to be found in paid subscriptions. I have defined a data_catalogue as a dictionary below which will contain symbols, data sources, and the like so our code knows where to go to get the data.

data_catalogue = {'DJIA':{
    'source': 'csv',
    'symbol': 'DJIA',
    'url': 'https://stooq.com/q/d/l/?s=^dji&i=d'
    },
    'S&P500': {
        'source': 'csv',
        'symbol': 'SPX',
        'url': 'https://stooq.com/q/d/l/?s=^spx&i=d'
    },
    'WTI': {
        'source': 'yahoo',
        'symbol': 'CL=F',
    },
    'Gold': {
        'source': 'yahoo',
        'symbol': 'GC=F',
    },
    'GBP': {
        'source': 'yahoo',
        'symbol': 'GBPUSD=X'
    },
    'BTC': {
        'source': 'yahoo',
        'symbol': 'BTC-USD'
    }
}

Now we’ll tie all of this together into a TestBench class. This will take our data catalogue, reshape it, and run our tests. The results are going to be collected for analysis, and I wrote a helper function to organize it into a large, Pandas dataframe for easy viewing.

import pandas as pd
import pandas_datareader as pdr
import yfinance as yf
from datetime import datetime

class TestBench:

  data_catalogue = data_catalogue

  test_names = ['runs-test',
                'dft-test',
                'bmr-test']

  def __init__(self, p_threshold=0.05, seed=101, 
               dftThreshold=0.95, bmrRows=32):
    np.random.seed(seed)
    self.seed = seed
    self.p_threshold = p_threshold
    self.dftThreshold = dftThreshold
    self.bmrRows = bmrRows
    self.years = [1, 4, 7, 10]
    self.trading_days = 250
    self.instruments = list(self.data_catalogue.keys())

  def getData(self):
    self.data_dict = {}
    for instr in self.instruments:
      try:
        data = self._getData(instr)
      except Exception as e:
        print(f'Unable to load data for {instr}')
        continue
      self.data_dict[instr] = data.copy()
    self.data_dict['baseline'] = np.random.binomial(1, 0.5, 
      size=self.trading_days * max(self.years) * 10)

  def _getData(self, instr):
    source = self.data_catalogue[instr]['source']
    sym = self.data_catalogue[instr]['symbol']
    if source == 'yahoo':
      return self._getYFData(sym)
    elif source == 'csv':
      return self._getCSVData(self.data_catalogue[instr]['url'])
    elif source == 'fred':
      return self._getFREDData(sym)

  def _getCSVData(self, url):
    data = pd.read_csv(url)
    close_idx = [i 
      for i, j in enumerate(data.columns) if j.lower() == 'close']
    assert len(close_idx) == 1, f"Can't match column names.\n{data.columns}"
    try:
      std_data = self._standardizeData(data.iloc[:, close_idx[0]])
    except Exception as e:
      raise ValueError(f"{url}")
    return std_data

  def _getYFData(self, sym):
    yfObj = yf.Ticker(sym)
    data = yfObj.history(period='max')
    std_data = self._standardizeData(data)
    return std_data


  def _getFREDData(self, sym):
    data = pdr.DataReader(sym, 'fred')
    data.columns = ['Close']
    std_data = self._standardizeData(data)
    return std_data


  def _standardizeData(self, df):
    # Converts data from different sources into np.array of price changes
    try:
      return df['Close'].diff().dropna().values
    except KeyError:
      return df.diff().dropna().values

  def runTests(self):
    self.test_results = {}
    for k, v in self.data_dict.items():
      print(f"{k}")
      self.test_results[k] = {}
      for t in self.years:
        self.test_results[k][t] = {}
        data = self._reshapeData(v, t)
        if data is None:
          # Insufficient data
          continue

        self.test_results[k][t]['runs-test'] = np.array(
          [self._runsTest(x) for x in data])
        self.test_results[k][t]['dft-test'] = np.array(
          [self._dftTest(x) for x in data])
        self.test_results[k][t]['bmr-test'] = np.array(
          [self._bmrTest(x) for x in data])

        print(f"Years = {t}\tSamples = {data.shape[0]}")

  def _reshapeData(self, X, years):
    d = int(self.trading_days * years) # Days per sample
    N = int(np.floor(X.shape[0] / d)) # Number of samples
    if N == 0:
      return None
    return X[-N*d:].reshape(N, -1)

  def _dftTest(self, data):
    return DFTTest(data, self.dftThreshold)

  def _runsTest(self, data):
    return RunsTest(data)

  def _bmrTest(self, data):
    return binMatrixRankTest(data, self.bmrRows)

  def tabulateResults(self):
    # Tabulate results
    table = pd.DataFrame()
    row = {}
    for k, v in self.test_results.items():
      row['Instrument'] = k
      for k1, v1 in v.items():
        row['Years'] = k1
        for k2, v2 in v1.items():
          pass_rate = sum(v2>self.p_threshold) / len(v2) * 100
          row['Test'] = k2
          row['Number of Samples'] = len(v2)
          row['Pass Rate'] = pass_rate
          row['Mean P-Value'] = v2.mean()
          row['Median P-Value'] = np.median(v2)
          table = pd.concat([table, pd.DataFrame(row, index=[0])])
    return table

We can initialize our test bench and call the getData() and runTests() method to put it all together. The tabulateResults() method will give us a nice table for viewing.

When we run our tests, we have a print out for the number of years and full samples of data we have. You’ll notice that for some of these (e.g. Bitcoin) we just don’t have a great amount of data to go off of, but we’ll do our best with what we do have.

tests = TestBench()
tests.getData()
tests.runTests()
DJIA
Years = 1	Samples = 129
Years = 4	Samples = 32
Years = 7	Samples = 18
Years = 10	Samples = 12
S&P500
Years = 1	Samples = 154
Years = 4	Samples = 38
Years = 7	Samples = 22
Years = 10	Samples = 15
WTI
Years = 1	Samples = 21
Years = 4	Samples = 5
Years = 7	Samples = 3
Years = 10	Samples = 2
Gold
Years = 1	Samples = 20
Years = 4	Samples = 5
Years = 7	Samples = 2
Years = 10	Samples = 2
GBP
Years = 1	Samples = 18
Years = 4	Samples = 4
Years = 7	Samples = 2
Years = 10	Samples = 1
BTC
Years = 1	Samples = 10
Years = 4	Samples = 2
Years = 7	Samples = 1
Years = 10	Samples = 1
baseline
Years = 1	Samples = 100
Years = 4	Samples = 25
Years = 7	Samples = 14
Years = 10	Samples = 10

We have 129 years of Dow Jones data, which gives us 12, 10-year samples and 154 years for the S&P 500 (the index doesn’t go back that far, but our data source provides monthly data going back to 1789). This is in contrast to most of our other values which have two decades or less.

To take a look at the results, we can run the tabulateResults() method, and do some pivoting to reshape the data frame for easier viewing.

table = tests.tabulateResults()
pivot = table.pivot_table(index=['Instrument', 'Years'], columns='Test')
samps = pivot['Number of Samples'].drop(['bmr-test', 'dft-test'], axis=1)
pivot.drop(['Number of Samples'], axis=1, inplace=True)
pivot['Number of Samples'] = samps
pivot

Let’s start with the baseline.

As expected, NumPy’s random number generator is pretty good, and it passes most of the tests without issue. The median P-values for the runs and DFT tests remain fairly high as well, although they are lower for the BMR test. Another thing to note, the 1 and 4 year BMR tests didn’t return any values because we were unable to complete a single 32×32 matrix with such small sample sizes. Overall, the lack of data for the BMR test makes the results here dubious (we could recalculate it with a smaller matrix size, but we’d need to recalibrate all of the probabilities for these different matrices).

The DFT test showed randomness for most cases in our test set. For what it’s worth, the P-values for our DFT tests of all sizes remained fairly high regardless of the sample size.

The runs test provides the most varied and interesting results.

import matplotlib.pyplot as plt

plt.figure(figsize=(12, 8))
for i, instr in enumerate(tests.instruments):
  sub = table.loc[(table['Instrument']==instr) &
                  (table['Test']=='runs-test')]
  plt.plot(tests.years, sub['Pass Rate'], label=instr, 
           c=colors[i], marker='o')

plt.legend()
plt.xlabel('Years')
plt.ylabel('Pass Rate (%)')
plt.title('Runs Test Pass Rate for all Instruments')
plt.show()

The runs test tends to produce less random results as time goes on. The notable exception being our WTI data, which passes more tests for randomness over time. However, if we look at our P-values, we do see them falling towards 0 (recall, our null hypothesis is that these are random processes).

plt.figure(figsize=(12, 8))
for i, instr in enumerate(table['Instrument'].unique()):
  sub = table.loc[(table['Instrument']==instr) &
                  (table['Test']=='runs-test')]
  plt.plot(tests.years, sub['Median P-Value'], label=instr, 
           c=colors[i], marker='o')
  
plt.legend()
plt.xlabel('Years')
plt.ylabel('P-Value')
plt.title('Median P-Values for Runs Test for all Instruments')
plt.show()

We added the baseline to this plot to show that it remains high even as the time frame increases, whereas all other values become less random over time. We’re showing P-values here, which are the probabilities that the results are due to noise if the process we’re testing is random. In other words, the lower our values become, the less likely it is that we have a random process on our hands.

This downward sloping trend may provide evidence that supports the value of longer-term trading.

Jerry Parker, for example, has moved toward longer-term trend signals (e.g. >200 day breakouts) because the short term signals are no longer profitable in his system. Data is going to be limited, but it could be interesting to run this over multiple, overlapping samples as in a walk forward analysis to see if randomness in the past was lower during shorter time frames. Additionally, there are more statistical tests we could look at to try to tease this out.

Death of the Random Walk Hypothesis?

The evidence from these few tests is mixed. Some tests show randomness, others provide an element of predictability. Unfortunately, we can’t definitively say the RWH is dead (although I think it, and the theories it is based on, are more articles of academic faith than anything).

To improve our experiment we need more data and more tests. We also used a series of binary tests, although technically the RWH asserts that the changes in price are normally distributed, so statistical tests that look for these patterns could strengthen our methodology and lead to more robust conclusions.

If you’d like to see more of this, drop us a note at hello@raposa.trade and let us know what you think!

Why I Left Value Investing Behind

You have probably heard the old adage “buy low and sell high.”

That’s great, but the question is what is “low” and what is “high?”

To know the difference between high and low, you need an evaluation framework. Here we have two main camps — the fundamental analysis camp (also known as value investing methodology) and the quantitative analysis camp.

Both methodologies are used to grow portfolios. I started off in the fundamental value camp in my early 20’s, studying the great value investors. I sorted through scores of annual reports, balance sheets, and income statements to build investment theses. Things worked out well, blessed by excellent timing, I was more lucky than good. No matter how many reports I read or industries I analyzed, I couldn’t separate the emotional and subjective nature of this approach from my decision making. On top of that, value has continued to lag other approaches since the 2008 crash.

Over time, I began building algorithms to guide my trading decisions. The more progress I made the more I transitioned into the quantitative analysis camp. The result of this journey has become Raposa Technologies — a way to validate your trading strategies without any coding knowledge.

Fundamental Analysis: One Stock at a Time

Fundamental analysis seeks to find the “intrinsic value” of a security based on attributes like free cash flow, debt levels, book value, and so forth. To do it well — like Warren Buffett — you need to read a copious amount of annual reports, quarterly earnings, and understand drivers in an industry you’re investing in. You will become an expert about companies you consider investing in. A core step of fundamental analysis is estimating all future cash flows and year over year profitability. In order to do this well— you must discount all future cash flows because money today is worth more than money 10 years from now.

So if you have the time to apply careful analysis to dozens of companies, read hundreds of reports, understand several industries, and carefully calculate future cash flows — — you will have an estimate for the fundamental price of a security. If the market price is lower than what your fundamental analysis estimates it to be, congratulations! You now have a good candidate to buy low while you wait for the market to realize the value of this stock and raise the price.

Value investors tend to get a lot of press (who hasn’t heard of Buffett?) because they can weave a narrative around a stock’s price journey. These narratives appeal to the emotional centers in our brains. Our brains are quick to use these stories as rationalization to support our gut feeling telling us to buy (or sell) a particular stock.

Your gut is quickly convinced by fundamental value narratives — particularly when they come from people who made fortunes riding these stocks to the top. Stories of double or triple digit returns from Amazon, Apple, Google, and even meme-stocks make it all too easy to believe the first fundamental narrative we hear.

During a bull market — it is easy to imagine the profits rolling in — but do not forget the emotional toll of holding names like Amazon through the long dark periods of doubt and uncertainty. You forget that their triumph wasn’t inevitable in the mid-2000s when the competitive landscape was forming. Could you have held on through the tech bust? What about the 2008 crash? Will you be confident in your fundamental analysis the morning you wake up to a 50–80% drop in your portfolio?

PROBABLY NOT.

But that’s fine. It takes a LOT of research and emotional work to invest in stocks based on fundamental analysis — which is why Buffett himself recommends people just buy an index fund and let it ride.

Quantitative Analysis: Separate the Signal from the Noise

After years of trying to invest based on fundamentals, not only was I treading water trying to balance my time — but I came to the realization that most of my investment decisions boiled down to a gut feeling no matter how rational and logical I tried to be.

There are successful quants and successful value investors. But if you are reading this far, you are probably unimpressed with the prospect of spending hundreds of hours researching companies for your fundamental analysis spreadsheets. You want new strategies in your war chest.

Quants are unconcerned about the intrinsic value of a stock or security, we look at the statistical profile of its price. How is it correlated with other prices? How does volume impact it? Are there regular patterns in price that can be leveraged for profit? Once we find a pattern — we can design algorithms to automatically execute trades that over time will grow our profiles.

These patterns often make no sense from a value perspective. Why would you buy a stock that appears to be incredibly overvalued? If you’re running a momentum or trend following strategy, you could find yourself buying near all-time highs. The value investor views that as insanity, but you do it because the algorithm shows that you have a potentially profitable pattern in your data set. That means you’re playing the odds that you can buy high and sell higher.

Do most investors have data-backed confidence for their trades? Or are decisions the results of a gut feeling? Considering most people run from math and code, I wager many trades are emotionally driven.

Break Away from the Narrative

Quantitative methods involve complicated statistical analysis, calculus, and machine learning capabilities. If you want to do it yourself, you’re going to need to learn to code. The upside? Your algorithms are working for you — which won’t eliminate your emotions or temptation to intervene — but the emotionless data will provide a beacon of rationality when FOMO or headline panic sets in.

For me, this was a big upside. I decided to apply my data science skills (those same skills I had honed during a PhD and applied everyday in a 9–5 for many years) and found that the math and stats in the quant world were much better for me, and improved my returns.

I firmly planted myself in this camp and never looked back.

I realize too that these methods aren’t easy, so that’s why I started Raposa — to build a no-code platform to enable investors to test and trade a variety of quantitative strategies. If you hate math and stats, then it’s not for you. Otherwise, join the waitlist and sign up below.

How To Reduce Lag In A Moving Average

Moving average indicators are commonly used to give traders a general idea about the direction of the trend by smoothing the price series. One of the big drawbacks to most common moving averages is the lag with which they operate. A strong trend up or down may take a long time to get confirmation from the series leading to lost profit.

In 2005, Alan Hull devised the Hull Moving Average (HMA) to address this problem.

The calculation is relatively straightforward and can be done in 4 steps after choosing the number of periods, N, to use in the calculation:

  1. Calculate the simple moving average over the past N periods.
    • SMA1 = SMA(price, N)
  2. Calculate the simple moving average over the past N/2 periods, rounded to the nearest whole value.
    • SMA2 = SMA(price, int(N/2))
  3. Multiply the shorter moving average by 2 and then subtract the first moving average from this.
    • SMA_diff = 2 * SMA2 - SMA1
  4. Take the moving average of this value over a period length equal to the square root of N, rounded to the nearest whole number.
    • HMA = SMA(SMA_diff, int(sqrt(N)))

This winds up being more responsive to recent changes in price because we’re taking the most recent half of our data and multiplying it by 2. This provides an additional weighting on those values before we smooth things out again with the final moving average calculation. Confusingly, many blogs list each of these moving averages as weighted moving averages, but never specify the weights themselves. Don’t worry about that, all we have are a few simple moving averages which are weighted before being combined at the end.

For completeness, we can also write this out mathematically.

If we are calculating the SMA at time t over the last N periods, we’re going to call this SMA^N_t. For moving averages, we’re just getting a summation over the last N prices (we’ll use P for prices) and dividing by N like so:

SMA_t^N = \frac{1}{N}\sum_{i=1}^{N} P_{i-N}

SMA_t^M = \frac{1}{M}\sum_{i=1}^{M} P_{i-M}

HMA_t^H = \frac{1}{H} \sum_{i=1}^{H} (2SMA^M_t - SMA^N_t)

where the symbols M and H are N/2 and the square root of N rounded to the nearest integer values.

M = \bigg\lfloor \frac{N}{2} \bigg\rceil

H = \bigg\lfloor \sqrt{N} \bigg\rceil

Hopefully, that’s all pretty straightforward. Let’s get to some examples in Python to illustrate how this works.

Hull Moving Average in Python

Like usual, let’s grab a few packages.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import yfinance as yf

From here, we can write a function to calculate the HMA in just three lines of code, corresponding to the three equations we showed above.

def calcHullMA(price: pd.Series, N=50):
  SMA1 = price.rolling(N).mean()
  SMA2 = price.rolling(int(N/2)).mean()
  return (2 * SMA2 - SMA1).rolling(int(np.sqrt(N))).mean()

We have our two moving averages, take the difference, and then smooth out the results with a third moving average. This function assumes we’re working with a Pandas data series and takes advantage of many of the methods that enables. Just be careful not to pass it a list or a NumPy array!

Getting Some Data

Let’s illustrate how this works on some historical data. I’m just getting a year’s worth from a common stock, DAL.

ticker = 'DAL'
start = '2014-01-01'
end = '2015-01-01'
yfObj = yf.Ticker(ticker)
data = yfObj.history(start=start, end=end)
data.drop(['Open', 'High', 'Low', 'Volume', 'Dividends', 'Stock Splits'],
          axis=1, inplace=True)

# Applying our function
N = 50
data[f'HMA_{N}'] = calcHullMA(data['Close'], N)

Take a look to see how it behaves:

plt.figure(figsize=(12, 8))
plt.plot(data['Close'], label='Close')
plt.plot(data[f'HMA_{N}'], label='HMA')
plt.xlabel('Date')
plt.ylabel('Price ($)')
plt.title(f'HMA and Price for {ticker}')
plt.legend()
plt.show()

As you can see, the HMA follows pretty closely. Of course, there is a lag as can be seen with some of the larger peaks and valleys over this time frame. Does it smooth well and with a lower lag than other moving averages as Hull intends?

To find out, let’s compare it to a typical, simple moving average and an exponential moving average (EMA). Like the HMA, the EMA is designed to be more responsive to recent price changes.

The code for the EMA calculation below was taken from a previous post you can dive into for further details.

def _calcEMA(P, last_ema, N):
    return (P - last_ema) * (2 / (N + 1)) + last_ema
    
def calcEMA(data: pd.DataFrame, N: int, key: str = 'Close'):
    # Initialize series
    sma = data[key].rolling(N).mean()
    ema = np.zeros(len(data)) + np.nan
    for i, _row in enumerate(data.iterrows()):
        row = _row[1]
        if np.isnan(ema[i-1]):
            ema[i] = sma[i]
        else:
            ema[i] = _calcEMA(row[key], ema[i-1], N)
    return ema

Plotting the results:

data[f'EMA_{N}'] = calcEMA(data, N)
data[f'SMA_{N}'] = data['Close'].rolling(N).mean()

plt.figure(figsize=(12, 8))
plt.plot(data['Close'], label='Close', linewidth=0.5)
plt.plot(data[f'HMA_{N}'], label='HMA')
plt.plot(data[f'EMA_{N}'], label='EMA')
plt.plot(data[f'SMA_{N}'], label='SMA')
plt.xlabel('Date')
plt.ylabel('Price ($)')
plt.title('Comparing 50-Day Moving Averages to Price')
plt.legend()
plt.show()

The plot looks pretty good. The HMA seems to track the price more closely than the other indicators while providing some good smoothing. However, we aren’t technical traders here at Raposa, so we need to do more than just look at a chart. We want to see the data!

To get an idea for the tracking error, we’re going to use the root mean square error (RMSE) to measure the difference between the indicator value and the price.

The RMSE is a common error metric that punishes deviations by squaring the error term. This means an error of 2 is 4 times greater than an error of 1! These squared errors all get summed up and then we take the square root of the values divided by the number of observations, n.

RMSE = \sqrt{\frac{\sum_t \big(\hat{P}_t - P_t \big)^2}{n}}

We’ll run our errors through a quick RMSE function we’ll write and see the results.

# Calculate tracking error
def calcRMSE(price, indicator):
  sq_error = np.power(indicator - price, 2).sum()
  n = len(indicator.dropna())
  return np.sqrt(sq_error / n)

hma_error = calcRMSE(data['Close'], data[f'HMA_{N}'])
ema_error = calcRMSE(data['Close'], data[f'EMA_{N}'])
sma_error = calcRMSE(data['Close'], data[f'SMA_{N}'])

print('Lag Error')
print(f'\tHMA = \t{hma_error:.2f}')
print(f'\tEMA = \t{ema_error:.2f}')
print(f'\tSMA = \t{sma_error:.2f}')
Lag Error
	HMA = 	1.65
	EMA = 	1.24
	SMA = 	1.53

Whoa! The HMA actually has greater error vs the price it’s tracking than the EMA and the SMA. This seems to cut against the intent of the HMA.

This is a small sample size, however, so maybe it really does have less lag than the other indicators and we just chose a bad stock and/or time frame.

Let’s test this by calculating the RMSE all of the stocks in the S&P 500 over the course of a year. Additionally, we’ll do this for different values of N to see if there’s any relationship between shorter or longer term values and the error.

Below, we have a helper function to calculate these values for us.

def calcErrors(data: pd.DataFrame, N: list):
  hma_error, sma_error, ema_error = [], [], []
  for n in N:
    hma = calcHullMA(data['Close'], n)
    ema = pd.Series(calcEMA(data, n), index=data.index)
    sma = data['Close'].rolling(n).mean()
    hma_error.append(calcRMSE(data['Close'], hma))
    ema_error.append(calcRMSE(data['Close'], ema))
    sma_error.append(calcRMSE(data['Close'], sma))
  
  return hma_error, ema_error, sma_error

The calcErrors function takes our data and a list of time periods to calculate the HMA, EMA, and SMA. From there, we calculate the RMSE for each series versus our closing price and return lists of each.

Next, we’ll loop over all the stocks in the S&P 500 and get the data for each. We’ll pass this to our error calculation function and collect the errors for each symbol.

We’re relying on the list of stocks in Wikipedia, which doesn’t necessarily correspond to how the symbols are represented in yfinance (e.g. Berkshire Hathaway has two classes of shares A’s and B’s, which cause issues) so we need to wrap this in a try-except statement for those edge cases. We’ll still get enough that we should be able to get a decent estimate.

# Sample 10 tickers from S&P 500
url = 'https://en.wikipedia.org/wiki/List_of_S%26P_500_companies'
table = pd.read_html(url)
df = table[0]
syms = df['Symbol']
start = '2019-01-01'
end = '2020-01-01'
N = [5, 10, 15, 20, 30, 50, 100, 150, 200]
for i, s in enumerate(syms):
  try:
    yfObj = yf.Ticker(s)
    data = yfObj.history(start=start, end=end)
  except:
    continue
  he, ee, se = calcErrors(data, N)
  if i == 0:
    hma_error = np.array(he)
    ema_error = np.array(ee)
    sma_error = np.array(se)
  else:
    hma_error = np.vstack([hma_error, he])
    ema_error = np.vstack([ema_error, ee])
    sma_error = np.vstack([sma_error, se])

# Drop rows with missing values
hma_error = hma_error[~np.isnan(hma_error).any(axis=1)]
ema_error = ema_error[~np.isnan(ema_error).any(axis=1)]
sma_error = sma_error[~np.isnan(sma_error).any(axis=1)]

After a few minutes, we can take a look at the mean tracking error across all of our metrics and tickers below:

Here we see that the HMA does track the price much better than other moving average measurements. There’s much less difference in short-time frames, but the values do start to diverge from one another fairly quickly and become more pronounced over time.

Trading with the Hull Moving Average

We could be more rigorous by tracking the deviation of our error measurements and getting more data, however for most purposes, it does seem as if the HMA does deliver on its promise to reducing lag. How do you trade it though?

The nice thing about the HMA, is that you can use it anywhere you’d use a moving average of any variety. You could build a whole new strategy around it, or just plug it into an existing system to see if you get any boost in your results.

We make all of that as easy as possible at Raposa, where we’re building a platform to allow you to backtest your ideas in seconds, without writing a single line of code. You can check out our free demo here!

How to Trade the MACD: Four Strategies with Backtests

The Moving Average Convergence-Divergence (MACD) is a popular and versatile indicator that appears in a number of trading systems. In it’s most basic form, we have the difference between two exponential moving averages (EMA), one fast and the other slow. The MACD is the difference between these two EMAs. From this simple beginning a host of other indicators such as signal lines and MACD bars are built. We’ll show you how to implement each of these, but refer back to this article for a thorough explanation.

In this post, we’re going to implement and backtest three versions of the MACD in Python to illustrate how to trade it and discuss areas for improvement.

Of course, if you want to skip all the code, you can try our new platform for free here to test these ideas.

Mean-Reverting MACD

A good place to get started with the MACD is to see when the signal diverges, i.e. when the MACD moves far away from 0. A mean-reverting strategy can quickly be built on top of this.

Our strategy is going to use standard parameters for the MACD. Our fast EMA will look over the past 12 days, and the slow EMA will run over 26 days. The model is going to buy our stock when the price breaks below a certain threshold, and will sell when the MACD converges back to 0. If the MACD runs high, we’ll short the stock and sell when it gets back to 0. We’re simply trying to jump on large, outlier movements in the price with the hope that the price will move back towards the longer EMA.

To get going, fire up Python and import the following packages.

import numpy as np
import pandas as pd
import yfinance as yf
import matplotlib.pyplot as plt

To start, we need to calculate the MACD. As stated above, the MACD is built on top of the EMA, so we’re going to write a few functions to calculate the EMA, one to plug in the values, the other that will initialize it and apply it to our data frame. From there, we can write our MACD function that will take these parameters, calculate the EMAs, our MACD, and fill in our data frame.

The code for the following functions was taken from this, previous post on the MACD and will go through all of the details on the code and calculations.

def _calcEMA(P, last_ema, N):
  return (P - last_ema) * (2 / (N + 1)) + last_ema

def calcEMA(data: pd.DataFrame, N: int, key: str = 'Close'):
  # Initialize series
  data['SMA_' + str(N)] = data[key].rolling(N).mean()
  ema = np.zeros(len(data)) + np.nan
  for i, _row in enumerate(data.iterrows()):
    row = _row[1]
    if np.isnan(ema[i-1]):
      ema[i] = row['SMA_' + str(N)]
    else:
      ema[i] = _calcEMA(row[key], ema[i-1], N)
  
  data['EMA_' + str(N)] = ema.copy()
  return data

def calcMACD(data: pd.DataFrame, N_fast: int, N_slow: int):
  assert N_fast < N_slow, 
    ("Fast EMA must be less than slow EMA parameter.")
  # Add short term EMA
  data = calcEMA(data, N_fast)
  # Add long term EMA
  data = calcEMA(data, N_slow)
  # Subtract values to get MACD
  data['MACD'] = data[f'EMA_{N_fast}'] - data[f'EMA_{N_slow}']
  # Drop extra columns
  data.drop(['Open', 'High', 'Low', 'Volume',
    'Dividends', 'Stock Splits'], axis=1, inplace=True)
  return data

Now with the MACD in place, we’re going to write two more functions. One is going to be our strategy function, and the other will calculate all of our returns. We’re breaking this second one into its own function because we can re-use the code later when we write other MACD strategies.

The MACDReversionStrategy is where our strategy is implemented (if you didn’t guess that by the name). This is a simple, vectorized implementation that is just going to look for our level value, and trade if the MACD moves outside of its bounds.

def calcReturns(df):
  df['returns'] = df['Close'] / df['Close'].shift(1)
  df['log_returns'] = np.log(df['returns'])
  df['strat_returns'] = df['position'].shift(1) * df['returns']
  df['strat_log_returns'] = df['position'].shift(1) * \
    df['log_returns']
  df['cum_returns'] = np.exp(df['log_returns'].cumsum()) - 1  
  df['strat_cum_returns'] = np.exp(
    df['strat_log_returns'].cumsum()) - 1
  
  return df

def MACDReversionStrategy(data, N_fast=12, N_slow=26, 
  level=1, shorts=True):
  df = calcMACD(data, N_fast, N_slow)
  # Drop extra columns
  df.drop(['Open', 'High', 'Low', 'Volume',
    'Dividends', 'Stock Splits'], axis=1, inplace=True) 
  df['position'] = np.nan
  df['position'] = np.where(df['MACD']<level, 1, df['position'])
  if shorts:
    df['position'] = np.where(df['MACD']>level, -1, df['position'])
  
  df['position'] = np.where(df['MACD'].shift(1)/df['MACD']<0, 
    0, df['position'])
  df['position'] = df['position'].ffill().fillna(0)
  
  return calcReturns(df)

The next step is to get data. We’ll rely on the yfinance package to get some free data from Yahoo! Finance. For backtests, more data is better, so we’re going to grab an older company so we can see how this performs over a long, multi-decade time horizon, so let’s choose 3M (ticker: MMM) and see how this strategy works.

ticker = 'MMM'
start = '2000-01-01'
end = '2020-12-31'
yfObj = yf.Ticker(ticker)
df = yfObj.history(start=start, end=end)
N_fast = 12
N_slow = 26
df_reversion = MACDReversionStrategy(df.copy(), 
  N_fast=N_fast, N_slow=N_slow, level=1)
fig, ax = plt.subplots(2, figsize=(12, 8), sharex=True)
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
ax[0].plot(df_reversion['strat_cum_returns']*100,
  label='MACD Reversion')
ax[0].plot(df_reversion['cum_returns']*100, label=f'{ticker}')
ax[0].set_ylabel('Returns (%)')
ax[0].set_title(
  'Cumulative Returns for MACD Reversion Strategy' +
  f'and Buy and Hold for {ticker}')
ax[0].legend()
ax[1].plot(df_reversion['MACD'])
ax[1].axhline(level, label='Short Level', color=colors[1],
  linestyle='--')
ax[1].axhline(-level, label='Long Level', color=colors[1],
  linestyle='--')
ax[1].axhline(0, label='Neutral', color='k', linestyle='--')
ax[1].set_xlabel('Date')
ax[1].set_ylabel('MACD')
ax[1].set_title(f'{N_fast}/{N_slow} MACD for {ticker}')
plt.tight_layout()
plt.show()

Our MACD mean reversion strategy outpaces the underlying stock from 2000–2016, before the strategy begins giving back some value and then underperforms for a few years. Starting in 2018 and continuing through the end of 2020, however, the strategy simply takes off!

One thing to notice about this model, is that the range of values for the MACD increases over time. It’s likely the case that the signal value we chose (-1 and 1 for simplicity) aren’t optimal and an adaptive value may prove better. The drawback of this is that we now have another parameter we need to fit, which can lead to overfitting our results.

MACD with Signal Line

This next method relies on the MACD signal line. This is just the EMA of the MACD itself. Because of this, we can re-use a lot of our code from the previous functions to add the signal line to our data.

def calcMACDSignal(data: pd.DataFrame, N_fast: int, N_slow: int,
  N_sl: int = 9):
  data = calcMACD(data, N_fast=N_fast, N_slow=N_slow)
  data = calcEMA(data, N_sl, key='MACD')
  # Rename columns
  data.rename(
    columns={f'SMA_{N_sl}': f'SMA_MACD_{N_sl}',
      f'EMA_{N_sl}': f'SignalLine_{N_sl}'}, inplace=True)
  return data

A common way to trade the MACD with the signal line, is to buy when the MACD is above the signal line, and short or sell if the lines cross again.

def MACDSignalStrategy(data, N_fast=12, N_slow=26, N_sl=9,
  shorts=False):
  df = calcMACDSignal(data, N_fast, N_slow, N_sl)
  df['position'] = np.nan
  df['position'] = np.where(df['MACD']>df[f'SignalLine_{N_sl}'], 1,
    df['position'])
  
  if shorts:
    df['position'] = np.where(df['MACD']<df[f'SignalLine_{N_sl}'],
      -1, df['position'])
  else:
    df['position'] = np.where(df['MACD']<df[f'SignalLine_{N_sl}'],   
      0, df['position'])
  df['position'] = df['position'].ffill().fillna(0)
  return calcReturns(df)

For the signal line, it’s fairly typical to look at the 9-day EMA, so that’s what we’ll use here. Again, no fiddling with the settings or optimization, we’re just taking a first pass at the backtest.

N_sl = 9
df_signal = MACDSignalStrategy(df.copy(), N_fast=N_fast,
  N_slow=N_slow, N_sl=N_sl)

fig, ax = plt.subplots(2, figsize=(12, 8), sharex=True)
ax[0].plot(df_signal['strat_cum_returns']*100,
  label='MACD Signal Strategy')
ax[0].plot(df_signal['cum_returns']*100, label=f'{ticker}')
ax[0].set_ylabel('Returns (%)')
ax[0].set_title(
  f'Cumulative Returns for MACD Signal Strategy and' + 
  f'Buy and Hold for {ticker}')
ax[0].legend()

ax[1].plot(df_signal['MACD'], label='MACD')
ax[1].plot(df_signal[f'SignalLine_{N_sl}'], 
  label='Signal Line', linestyle=':')
ax[1].set_xlabel('Date')
ax[1].set_ylabel('MACD')
ax[1].set_title(f'{N_fast}/{N_slow}/{N_sl} MACD for {ticker}')
ax[1].legend()

plt.tight_layout()
plt.show()

Here, our results weren’t nearly as great compared to the mean reversion model (the interplay between the signal line and the MACD is hard to discern at this time scale, zooming in will show its many crosses). This strategy does, however, avoid large or long drawdowns. It seems to be slow and steady, but does fail to win the race via a long-only strategy.

In a more complex strategy where you are managing a portfolio of strategies, an equity curve like we’re showing for the MACD signal strategy could be very valuable. The steady returns could provide a safe haven to allocate capital to if you’re high-return strategies suddenly underperform due to a shift in the market regime.

If you’re looking for all out returns, changing up some of the parameters could be a good way to go beyond the standard values we used here.

MACD Momentum

Another technique that is frequently used is to take the difference between the MACD and the signal line. This is then plotted as vertical bars and called MACD bar or MACD histogram plots. If the bars grow over time, then we have increasing momentum and go long, if they decrease, then we have slowing momentum and a possible reversal.

Calculating this is rather straightforward, we just take our MACD and subtract off the signal line. Again, we can re-use a lot of our previous code to keep things short and sweet.

def calcMACDBars(data: pd.DataFrame, N_fast: int, N_slow: int,
  N_sl: int = 9):
  data = calcMACDSignal(data, N_fast=N_fast, N_slow=N_slow,
    N_sl=N_sl)
  data['MACD_bars'] = data['MACD'] - data[f'SignalLine_{N_sl}']
  return data

The histogram is positive when the MACD is above the signal line, and negative when it falls below.

This is an indicator of an indicator and is a few steps removed from the actual price. Regardless, there are a variety of ways to trade using it. One way is to look for increasing or decreasing momentum. If consecutive, positive bars grow, it means that we have a bullish signal for increasing momentum as the MACD moves away from the signal line. If they’re moving closer to one another, then we have a bearish signal and can short it.

We can also look for crossovers, however this is the same as the signal line strategy we implemented above.

You’ll also see things such as peak-trough divergences or slant divergences to generate signals. These look for consecutive hills in the MACD histogram chart and are usually picked up via visual inspection.

We’re going to start with a momentum strategy that looks at the growth of consecutive bars. We’ll use the standard, 12/26/9 format for the MACD signal line parameters, but we’ll see if we can pick up on three consecutive days of growth (positive or negative) in the bars for buy/short signals.

def MACDMomentumStrategy(data, N_fast=12, N_slow=26, N_sl=9,
  N_mom=3, shorts=False):
  df = calcMACDBars(data, N_fast=N_fast, N_slow=N_slow, N_sl=N_sl)
  df['growth'] = np.sign(df['MACD_bars'].diff(1))
  df['momentum'] = df['growth'].rolling(N_mom).sum()
  if shorts:
    df['position'] = df['momentum'].map(
      lambda x: np.sign(x) * 1 if np.abs(x) == N_mom else 0)
  else:
    df['position'] = df['momentum'].map(
      lambda x: 1 if x == N_mom else 0)
  df['position'] = df['position'].ffill().fillna(0)
  return calcReturns(df)

Running and plotting:

N_mom = 3

df_mom = MACDMomentumStrategy(df.copy(), N_fast=N_fast, 
  N_slow=N_slow,N_sl=N_sl, N_mom=N_mom)

fig, ax = plt.subplots(3, figsize=(12, 8), sharex=True)
ax[0].plot(df_mom['strat_cum_returns']*100, label='MACD Momentum')
ax[0].plot(df_mom['cum_returns']*100, label=f'{ticker}')
ax[0].set_ylabel('Returns (%)')
ax[0].set_title(f'Cumulative Returns for MACD Signal Strategy and Buy and Hold for {ticker}')
ax[0].legend()
ax[1].plot(df_mom['MACD'], label='MACD')
ax[1].plot(df_mom[f'SignalLine_{N_sl}'], label='Signal Line', 
  linestyle=':')
ax[1].set_title(f'{N_fast}/{N_slow}/{N_sl} MACD for {ticker}')
ax[1].legend()
ax[2].bar(df_mom.index, df_mom['MACD_bars'], label='MACD Bars',
  color=colors[4])
ax[2].set_xlabel('Date')
ax[2].set_title(f'MACD Bars for {ticker}')
plt.tight_layout()
plt.show()

This one doesn’t look so great. 20 years and losing 15% is not a strategy I want to be invested in. If we look a little deeper into this one, we can see that we jump in and out of a lot of positions very quickly. Usually we have some momentum, but because the momentum signal we’re looking at can change after a day or two, we aren’t in a position long enough to ride a trend. Essentially, we’re acting like a trend follower, but we bail so quickly we wind up taking a lot of small losses and have very few gains to show for it.

To alleviate this, let’s see if we can improve by using our MACD bars to signal momentum and only exit a position when the MACD bar crosses over. To see this, we’ll look for a sign change from positive to negative or vice versa.

def MACDSignalMomentumStrategy(data, N_fast=12, N_slow=26, N_sl=9,
  N_mom=3, shorts=False):
  df = calcMACDBars(data, N_fast=N_fast, N_slow=N_slow, N_sl=N_sl)
  df['growth'] = np.sign(df['MACD_bars'].diff(1))
  df['momentum'] = df['growth'].rolling(N_mom).sum()
  # Enter a long/short position if momentum is going in the right
  # direction and wait for cross-over
  position = np.zeros(len(df)) + np.nan
  for i, _row in enumerate(data.iterrows()):
    row = _row[1]
    mom = row['momentum']
    if np.isnan(mom):
      last_row = row.copy()
      continue
    if np.abs(mom) == N_mom and position[i-1] == 0:
      # Enter new position
      if shorts:
        position[i] = np.sign(mom)
      else:
        position[i] = 1 if np.sign(mom) == 1 else 0
    elif row['MACD_bars'] / last_row['MACD_bars'] < 0:
      # Change in sign indicates cross-over -> exit
      position[i] = 0
    else:
      # Hold position
      position[i] = position[i-1]
  
  df['position'] = position
  
  return calcReturns(df)

Running this combined strategy:

df_sig_mom = MACDSignalMomentumStrategy(df.copy(), N_fast=N_fast, 
  N_slow=N_slow,N_sl=N_sl, N_mom=N_mom)
fig, ax = plt.subplots(3, figsize=(12, 8), sharex=True)
ax[0].plot(df_sig_mom['strat_cum_returns']*100, 
  label='MACD Momentum')
ax[0].plot(df_sig_mom['cum_returns']*100, label=f'{ticker}')
ax[0].set_ylabel('Returns (%)')
ax[0].set_title(f'Cumulative Returns for MACD Signal Strategy' +
  f'and Buy and Hold for {ticker}')
ax[0].legend()
ax[1].plot(df_sig_mom['MACD'], label='MACD')
ax[1].plot(df_sig_mom[f'SignalLine_{N_sl}'], label='Signal Line', 
  linestyle=':')
ax[1].set_title(f'{N_fast}/{N_slow}/{N_sl} MACD for {ticker}')
ax[1].legend()
ax[2].bar(df_sig_mom.index, df_sig_mom['MACD_bars'], 
  label='MACD Bars', color=colors[4])
ax[2].set_xlabel('Date')
ax[2].set_title(f'MACD Bars for {ticker}')
plt.tight_layout()
plt.show()

This combined model does perform better, yielding us about 250% over this time frame. That is still about half of the long-only strategy that we’re trying to beat, but we could continue to tinker with these ideas to come up with something even better.

Before moving on, let’s take a look at some of the key metrics for each of these strategies.

Comparing Strategies

Below, we use have a helper function to get our strategy statistics.

def getStratStats(log_returns: pd.Series, 
  risk_free_rate: float = 0.02):
  stats = {}
  # Total Returns
  stats['tot_returns'] = np.exp(log_returns.sum()) - 1
  # Mean Annual Returns
  stats['annual_returns'] = np.exp(log_returns.mean() * 252) - 1
  
  # Annual Volatility
  stats['annual_volatility'] = log_returns.std() * np.sqrt(252)
  # Sharpe Ratio
  stats['sharpe_ratio'] = (
    (stats['annual_returns'] - risk_free_rate) 
    / stats['annual_volatility'])
  # Max Drawdown
  cum_returns = log_returns.cumsum()
  peak = cum_returns.cummax()
  drawdown = peak - cum_returns
  stats['max_drawdown'] = drawdown.max()
  # Max Drawdown Duration
  strat_dd = drawdown[drawdown==0]
  strat_dd_diff = strat_dd.index[1:] - strat_dd.index[:-1]
  strat_dd_days = strat_dd_diff.map(lambda x: x.days).values
  strat_dd_days = np.hstack([strat_dd_days,
    (drawdown.index[-1] - strat_dd.index[-1]).days])
  stats['max_drawdown_duration'] = strat_dd_days.max()
  return stats

Getting stats for each strategy:

stats_rev = getStratStats(df_reversion['strat_log_returns'])
stats_sig = getStratStats(df_signal['strat_log_returns'])
stats_mom = getStratStats(df_mom['strat_log_returns'])
stats_sig_mom = getStratStats(df_sig_mom['strat_log_returns'])
stats_base = getStratStats(df_reversion['log_returns'])
stats_dict = {'Mean Reversion': stats_rev,
  'Signal Line': stats_sig,
  'Momentum': stats_mom,
  'Signal-Momentum': stats_sig_mom,
  'Baseline': stats_base}
pd.DataFrame(stats_dict)

This baseline was pretty tough to beat. Over 9% annual returns despite some large drawdowns, however, the MACD mean reversion strategy beat it easily and with a better Sharpe Ratio to boot. The signal line and the momentum model both severely underperformed the baseline — with the MACD momentum model just being flat out atrocious. However, combining these two yielded better results, still below the baseline, but with a reasonable Sharpe Ratio and a lower drawdown (albeit longer than the baseline).

No Guarantees

Each of these models could be tweaked and improved, but they’re only a small taste of the possibilities available to trade using the MACD and its derivatives. No matter how good that mean reversion strategy looks, nothing here should be blindly implemented.

While we ran a fine backtest, there are issues. Dividends weren’t taken into account, transaction costs were ignored, slippage, only one stock was examined, and so forth. If you really want to trade algorithmically, you’re better off with a strategy that manages a diversified portfolio, on data that has been adjusted for dividends, with transaction costs, over a long time horizon.

At Raposa, we’re building an engine that will enable you to quickly test ideas like these on a portfolio of securities, with proper money management, a plethora of signals, and high quality data, all without a single line of code. We have a state-of-the-art system that will allow you to customize your model to suit your needs, and deploy it when you’re ready so you can trade live. If you want to move beyond simple backtests and strategies, then sign up below to be kept up to date on our latest developments.

How to Build a Better Trading System

Most trading strategies you find online are of questionable value. So if you’re going to put your money into it, we suggest you build and test it yourself.

Thankfully, building your own strategy doesn’t have to be a daunting task – we provide tools to do just that – and outline 6 steps you can follow to build your own strategy that you can have confidence in.

Data or Ideas?

Profitable strategies can come from a variety of sources. Most algorithmic traders break it into idea-first and data-first strategies.

Idea-first strategies starts with a hypothesis about what might be profitable. For example, you might wonder if some indicator can be applied to a given stock or currency or if there’s a way to profit on cognitive biases. So you go ahead, get the data and test your idea.

Data-first strategies start with the data and try to extract potentially profitable patterns from the data which then form your hypothesis. These can be found using black-box models (e.g. machine learning) or oddball patterns you observe in the market, even if you can’t explain why they work.

Gregory Zuckerman explains this type of thinking while asking an executive at the famous quant fund Renaissance Technologies about signals they’d trade. As long as they had the statistics to back it up, they’d trade such strange signals as “volume divided by price change three days earlier;” it doesn’t matter if they have a story about why that works or not, just go for it.

Which is better?

Well, that’s hard to say – traders have had a lot of success with both approaches. The data-first approach is often more mathematically challenging, but that doesn’t mean it’s going to be more profitable. I agree with systematic trader and author Rob Carver when he writes, “consistently profitable trading comes out of careful research, done by thoughtful and knowledgeable people, who seek to understand where their profits come from. The loss-making systematic trading I’ve seen has often been the result of haphazard data mining…That experience, combined with my preference for things I can trust and understand, means I favor the ideas-first method.”

Test a Simple Version of Your Strategy

Now you have some idea, so you need to go test it to see if it really holds up. I suggest starting simple, without a lot of fancy position sizing, use of stop losses, checking correlations, and so forth – unless these are key parts of your idea! Just try to test the simplest version you can to see if there’s some potential.

You’re not looking for a world-beating backtest at this point, you just want to know whether or not there’s some potential in the signal you’re trying to exploit. Does it seem to do better in a bear market or bull market? What about high volatility regimes vs low volatility?

If there’s an edge in some situations, you might have something you can build and work with!

Add the Bells and Whistles

Assuming you’ve got a trading signal that piques your interest, you can start adding in some of the key components that a live strategy is going to need. You’re going to want to work with position sizing and risk management to avoid blowing up, add stops/targets to get out of trades, and add filters or other signals to restrict entry into a trade during periods where your model will perform at its best.

Optimize your Parameters

A lot of traders get into trouble by over-fitting their strategies. They jump on a signal and keep tweaking parameters until they get a model with an astronomical return. Lured in by the promise of riches, they don’t realize that their model is incredibly fragile and doomed to failure.

File:Overfitting.svg - Wikimedia Commons
The green line is over-fit to the data. It’s going to show great stats, but looking more closely, it’s clear that it is going to struggle with new data (image from Wikipedia).

While there are no hard and fast rules to avoid over-fitting, it’s a good idea to limit the number of parameters in your model and the number of runs you try.

Too many parameters allow you to play with a lot of combinations to find that combination that is “just right” and looks great in a backtest, but doesn’t generalize to your trading account. Each set of parameters requires a new run, so if you find yourself running “My Retire-Next-Year Strategy #149285” then it’s safe to say you should give it a break and try a new idea.

Out-of-Sample Testing

Let’s say you have 20 years of historical data available – most novice traders are going to fit their strategy on all 20 years, then go and trade. A better approach is to split between test and training data, so you optimize your parameters on the first 15 years, then test the results on the last 5 years. The first 15 years are your in-sample data while the last 5 years comprise your out-of-sample data. This helps prevent over-fitting as discussed above because you should be able to see how much your strategy’s performance degrades during the out-of-sample test.

Some degradation is expected, so your in-sample test should have higher returns and better risk metrics than your out-of-sample test. But if in-sample is amazing and out-of-sample is horrendous, then you’ve probably over-fit your data and need to go back to the drawing board.

A lot more can be said about proper testing and optimization (we’ll go into details in future posts). One of the best techniques to use is walk forward optimization. This is where you optimize on a small, subset of in-sample data (e.g. 1 year) then run a test on the next subset of your data for out of sample testing. You can do this with a lot of different parameters and keep the top 30% every time and see what survives. This requires a lot of data and discipline, but is widely considered the gold-standard approach.

Sanity Check

Congrats if your strategy has made it this far! Look closely at those returns though, do you really believe they’re possible in practice?

Carver argues that a single-instrument strategy should produce a realistic Sharpe Ratio of 0.3. You can get higher Sharpe Ratios using proper money and risk management as well as a broader portfolio of instruments, but if all of that gets well above 1, then you’re probably traipsing into fantasy land.

If you want to hold onto this strategy, then it may be prudent to paper trade for a bit while gathering more data to see if it performs the way you expected.

Go Trade!

Now you should have a clearly defined, optimized, and well tested trading system that provides attractive yet reasonable returns. Looks like it’s time to trade it!

Some people build custom dashboards, spreadsheets, and all sorts of infrastructure to manage and trade their own portfolio. An easier solution would be to work with a system that will provide you automated alerts when your time to trade has come up. Even better would be one that enables you to also quickly test a variety of ideas in an easy, no-code framework like this one here.

If you’re interested in learning more about Raposa and our cloud-based trading products, you can also join our list below!

4 Ways to Trade the Trend Intensity Indicator

Determining the strength of a trend can provide a valuable edge to your trading strategy and help you determine when to go long and let it ride, or not. This is what the Trend Intensity Indicator (TII) was designed to do.

This indicator is as simple to interpret as more familiar values like the RSI. It’s scaled from 0-100 where higher numbers indicate a stronger upward trend, lower values a stronger downward trend, and values around the centerline (50) are neutral.

Calculating the TII takes place in 3 quick steps:

  1. Calculate the Simple Moving Average (SMA) of your closing prices (60-day SMA is typical, i.e. P=60).
    • `SMA[t] = mean(Price[t-P:t])`
  2. Determine how many time periods close above the SMA over P/2 periods.
    • `PositivePeriods[t] = count(SMA[t-P/2:t]>Price[t-P/t:t]`
  3. Divide the number of positive periods by P/2 and multiply by 100 to scale from 0-100.
    • `TII[t] = PositivePeriods[t] / (P/2) * 100`

TL;DR

We provide backtests and code for four distinct trading strategies using the TII. These tests are designed to give you a feel for the indicator and running backtests. If you want something faster and more complete, check out our free, no-code platform here.

Trend Intensity Indicator

Let’s get to coding this in Python:

import pandas as pd

def calcTrendIntensityIndex(data, P=60):
  data[f'SMA_{P}'] = data['Close'].rolling(P).mean()
  diff = data['Close'] - data[f'SMA_{P}']
  pos_count = diff.map(lambda x: 1 if x > 0 else 0).rolling(int(P/2)).sum()
  data['TII'] = 200 * (pos_count) / P
  return data

The function is quick and easy.

Let’s grab some data and look a bit closer before trading it.

import numpy as np
import matplotlib.pyplot as plt
import yfinance as yf

ticker = 'GDX'
start = '2011-01-01'
end = '2013-01-01'
yfObj = yf.Ticker(ticker)
data = yfObj.history(start=start, end=end)
# Drop unused columns
data.drop(
    ['Open', 'High', 'Low', 'Volume', 'Dividends', 'Stock Splits'], 
    axis=1, inplace=True)

P = [30, 60, 90]

colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

fig, ax = plt.subplots(2, figsize=(12, 8), sharex=True)
ax[0].plot(data['Close'], label='Close')

for i, p in enumerate(P):
  df = calcTrendIntensityIndex(data.copy(), p)
  ax[0].plot(df[f'SMA_{p}'], label=f'SMA({p})')
  ax[1].plot(df['TII'], label=f'TII({p})', c=colors[i+1])

ax[0].set_ylabel('Price ($)')
ax[0].set_title(f'Price and SMA Values for {ticker}')
ax[0].legend()

ax[1].set_xlabel('Date')
ax[1].set_ylabel('TII')
ax[1].set_title('Trend Intensity Index')
ax[1].legend()
plt.tight_layout()
plt.show()

Just like the RSI and other oscillators, the TII ranges from 0-100. It does so in a steadier fashion than many of these other indicators because it’s simply counting the number of days where the price is greater than the SMA over the past P/2 periods (30 days in this case). It only moves in discrete steps up and down over time.

Even though the TII is derived from the SMA crossover, it differs significantly from the standard indicator as shown below.

P = 60
df = calcTrendIntensityIndex(data.copy(), P)
vals = np.where(df[f'SMA_{P}'] > df['Close'], 1, 0)
df['SMA_cross'] = np.hstack([np.nan, vals[:-1] - vals[1:]])

fig, ax = plt.subplots(2, figsize=(12, 8), sharex=True)
ax[0].plot(df['Close'], label='Close')
ax[0].plot(df[f'SMA_{P}'], label=f'SMA({P})')
ax[0].scatter(df.loc[df['SMA_cross']==1].index, 
              df.loc[df['SMA_cross']==1][f'SMA_{P}'], 
              marker='^', c=colors[3],
              label='Cross Above', zorder=100, s=100)
ax[0].scatter(df.loc[df['SMA_cross']==-1].index, 
              df.loc[df['SMA_cross']==-1][f'SMA_{P}'], 
              marker='v', c=colors[4],
              label='Cross Below', zorder=100, s=100)

ax[0].set_ylabel('Price ($)')
ax[0].set_title(f'Price and {P}-Day SMA for {ticker}')
ax[0].legend()

ax[1].plot(df['TII'], label='TII')
ax[1].scatter(df.loc[df['SMA_cross']==1].index,
              df.loc[df['SMA_cross']==1]['TII'], 
              marker='^', c=colors[3],
              label='Cross Above', zorder=100, s=100)
ax[1].scatter(df.loc[df['SMA_cross']==-1].index, 
              df.loc[df['SMA_cross']==-1]['TII'], 
              marker='v', c=colors[4],
              label='Cross Below', zorder=100, s=100)
ax[1].set_xlabel('Date')
ax[1].set_ylabel('TII')
ax[1].set_title('Trend Intensity Index')
ax[1].legend()
plt.tight_layout()
plt.show()

In this plot, we have the 60-day SMA marked with orange triangles when the price moves above the SMA and blue triangles showing when it drops below the SMA. As can be seen in the lower plot, these occur all over the TII. So you can’t necessarily tell when a cross over occurs by looking at the TII. While it is a straightforward derivative of the SMA crossover, you’re going to have very different signals when trading it.

Let’s turn to a few examples of how we can use it.

Buy Low and Sell High

The first model we’ll run approaches the TII in a similar fashion as the RSI; we buy when the TII crosses below 20 and sell/short when it breaks above 80. We hold the positions until they reach the exit value.

The basic idea is that the price should be split roughly between days above and below the SMA. When we have a series that has deviated from the long-run average, we put a position on and wait for it to revert back to the mean.

def TIIMeanReversion(data, P=60, enter_long=20, exit_long=50, 
                     enter_short=80, exit_short=50, shorts=True):
  df = calcTrendIntensityIndex(data, P=P)
  df['position'] = np.nan
  df['position'] = np.where(df['TII']<=enter_long, 1, np.nan)
  _exit_long = df['TII'] - exit_long
  exit = _exit_long.shift(1) / _exit_long
  df['position'] = np.where(exit<0, 0, df['position'])

  if shorts:
    df['position'] = np.where(df['TII']>=enter_short, -1, 
                              df['position'])
    _exit_short = df['TII'] - exit_short
    df['position'] = np.where(exit<0, 0, df['position'])

  df['position'] = df['position'].ffill().fillna(0)
  return calcReturns(df)

# A few helper functions
def calcReturns(df):
  # Helper function to avoid repeating too much code
  df['returns'] = df['Close'] / df['Close'].shift(1)
  df['log_returns'] = np.log(df['returns'])
  df['strat_returns'] = df['position'].shift(1) * df['returns']
  df['strat_log_returns'] = df['position'].shift(1) * df['log_returns']
  df['cum_returns'] = np.exp(df['log_returns'].cumsum()) - 1
  df['strat_cum_returns'] = np.exp(df['strat_log_returns'].cumsum()) - 1
  df['peak'] = df['cum_returns'].cummax()
  df['strat_peak'] = df['strat_cum_returns'].cummax()
  return df

def getStratStats(log_returns: pd.Series,
  risk_free_rate: float = 0.02):
  stats = {}  # Total Returns
  stats['tot_returns'] = np.exp(log_returns.sum()) - 1  
  
  # Mean Annual Returns
  stats['annual_returns'] = np.exp(log_returns.mean() * 252) - 1  
  
  # Annual Volatility
  stats['annual_volatility'] = log_returns.std() * np.sqrt(252)
  
  # Sortino Ratio
  annualized_downside = log_returns.loc[log_returns<0].std() * \
    np.sqrt(252)
  stats['sortino_ratio'] = (stats['annual_returns'] - \
    risk_free_rate) / annualized_downside  
  
  # Sharpe Ratio
  stats['sharpe_ratio'] = (stats['annual_returns'] - \
    risk_free_rate) / stats['annual_volatility']  
  
  # Max Drawdown
  cum_returns = log_returns.cumsum() - 1
  peak = cum_returns.cummax()
  drawdown = peak - cum_returns
  max_idx = drawdown.argmax()
  stats['max_drawdown'] = 1 - np.exp(cum_returns[max_idx]) \
    / np.exp(peak[max_idx])
  
  # Max Drawdown Duration
  strat_dd = drawdown[drawdown==0]
  strat_dd_diff = strat_dd.index[1:] - strat_dd.index[:-1]
  strat_dd_days = strat_dd_diff.map(lambda x: x.days).values
  strat_dd_days = np.hstack([strat_dd_days,
    (drawdown.index[-1] - strat_dd.index[-1]).days])
  stats['max_drawdown_duration'] = strat_dd_days.max()
  return {k: np.round(v, 4) if type(v) == np.float_ else v
          for k, v in stats.items()}

I also added a few helper functions to get our returns and statistics for easy plotting and comparison.

Let’s get some extra data for our ticker to see if we can hit it big more than our gold mining ETF.

ticker = 'GDX'
start = '2000-01-01'
end = '2020-12-31'
yfObj = yf.Ticker(ticker)
data = yfObj.history(start=start, end=end)
# Drop unused columns
data.drop(
    ['Open', 'High', 'Low', 'Volume', 'Dividends', 'Stock Splits'], 
    axis=1, inplace=True)

P = 60
enter_long = 20
enter_short = 80
df_mr = TIIMeanReversion(data.copy(), P=P, enter_long=enter_long, 
                         enter_short=enter_short)

fig, ax = plt.subplots(3, figsize=(12, 8), sharex=True)
ax[0].plot(df_mr['Close'], label='Close')
ax[0].plot(df_mr[f'SMA_{P}'], label=f'SMA({P})')
ax[0].set_ylabel('Price ($)')
ax[0].set_title(f'Price and {P}-Day SMA for {ticker}')
ax[0].legend()

ax[1].plot(df_mr['TII'], label='TII')
ax[1].axhline(enter_short, c=colors[1], 
              label='Short Signal Line', linestyle=':')
ax[1].axhline(enter_long, c=colors[2], 
              label='Long Signal Line', linestyle=':')
ax[1].axhline(50, c=colors[3],
              label='Exit Line', linestyle=':')
ax[1].set_ylabel('TII')
ax[1].set_title('Trend Intensity Index')
ax[1].legend(bbox_to_anchor=[1, 0.65])

ax[2].plot(df_mr['strat_cum_returns']*100, label='Mean Reversion')
ax[2].plot(df_mr['cum_returns']*100, label='Buy and Hold')
ax[2].set_title('Cumulative Returns')
ax[2].set_ylabel('Returns (%)')
ax[2].set_xlabel('Date')
ax[2].legend()
plt.tight_layout()
plt.show()

df_stats = pd.DataFrame(getStratStats(df_mr['log_returns']),
                        index=['Buy and Hold'])
df_stats = pd.concat([df_stats,
                      pd.DataFrame(getStratStats(
                        df_mr['strat_log_returns']),
                        index=['Mean Reversion'])])
df_stats

The mean reversion model yields some reasonable returns above and beyond the buy and hold approach. The risk adjusted metrics aren’t great, but it is a good boost against the underlying.

Momentum Trading

Like RSI and other oscillators, we can interpret the TII as a momentum indicator, buying when it breaks above the centerline and selling/shorting when it crosses below it. For this, we’ll add a bit of a wrinkle such that we’ll hold while the TII goes above an upper level and then sell if it proceeds to fall below that level. We treat the short side the same way. This will get us out of positions sooner, hopefully with a larger profit than if we held until it reversed all the way back to the centerline again.

def TIIMomentum(data, P=60, centerline=50, upper=80,
                lower=20, shorts=True):
  df = calcTrendIntensityIndex(data, P=P)
  position = np.zeros(df.shape[0])
  for i, (idx, row) in enumerate(df.iterrows()):
    if np.isnan(row['TII']):
      last_row = row.copy()
      continue
    
    if row['TII'] >= centerline and last_row['TII'] < centerline:
      # Go long if no position
      if position[i-1] != 1:
        position[i] = 1
      
    elif row['TII'] > centerline and position[i-1] == 1:
      # Check if broke below upper line
      if last_row['TII'] > upper and row['TII'] <= upper:
        position[i] = 0
      else:
        position[i] = 1
    
    elif position[i-1] == 1 and row['TII'] <= centerline:
      # Sell/short if broke below centerline
      if shorts:
        position[i] = -1
      else:
        position[i] = 0
    
    elif shorts:
      if row['TII'] <= centerline and last_row['TII'] > centerline:
        # Go short if no position
        if position[i-1] != -1:
          position[i] = -1
        elif row['TII'] <= centerline and position[i-1] == -1:
          # Check if broke above lower line
          if last_row['TII'] < lower and row['TII'] > lower:
            position[i] = 0
          else:
            position[i] = -1
        elif position[i-1] == -1 and row['TII'] > centerline:
          # Exit short and go long if it crosses
          position[i] = 1

    last_row = row.copy()
    
  df['position'] = position  
  
  return calcReturns(df)

The easiest way to get the logic right is just by looping over the data and buying/selling in each of the different cases. The code is a bit longer, but hopefully easier to understand.

Let’s move on to testing it with our standard settings.

P = 60
df_mom = TIIMomentum(data.copy(), P=P)

fig, ax = plt.subplots(3, figsize=(12, 8), sharex=True)
ax[0].plot(df_mom['Close'], label='Close')
ax[0].plot(df_mom[f'SMA_{P}'], label=f'SMA({P})')
ax[0].set_ylabel('Price ($)')
ax[0].set_title(f'Price and {P}-Day SMA for {ticker}')
ax[0].legend()

ax[1].plot(df_mom['TII'], label='TII')
ax[1].axhline(20, c=colors[1], 
              label='Exit Short Signal Line', linestyle=':')
ax[1].axhline(80, c=colors[2], 
              label='Exit Long Signal Line', linestyle=':')
ax[1].axhline(50, c=colors[3],
              label='Center Line', linestyle=':')
ax[1].set_ylabel('TII')
ax[1].set_title('Trend Intensity Index')
ax[1].legend(bbox_to_anchor=[1, 0.65])

ax[2].plot(df_mom['strat_cum_returns']*100, label='Momentum')
ax[2].plot(df_mom['cum_returns']*100, label='Buy and Hold')
ax[2].set_title('Cumulative Returns')
ax[2].set_ylabel('Returns (%)')
ax[2].set_xlabel('Date')
ax[2].legend()
plt.tight_layout()
plt.show()

df_stats = pd.concat([df_stats,
                      pd.DataFrame(getStratStats(
                        df_mom['strat_log_returns']),
                        index=['Momentum'])])
df_stats

This strategy lost over 60% of its initial capital – clearly not what we were hoping for. The underlying is fairly volatile, so it may perform better if we reduce P to try to take advantage of some of these shorter moves. I leave that as an exercise to the reader.

Adding a Signal Line

Like the MACD, the TII is often traded with a signal line. This signal line is just the EMA of the TII (i.e. EMA(TII(P))).
To calculate this, we need to grab a couple of functions for the EMA calculation.

def _calcEMA(P, last_ema, N):
  return (P - last_ema) * (2 / (N + 1)) + last_ema

def calcEMA(series, N):
  series_mean = series.rolling(N).mean()
  ema = np.zeros(len(series))
  for i, (_, val) in enumerate(series_mean.iteritems()):
    if np.isnan(ema[i-1]):
      ema[i] += val
    else:
      ema[i] += _calcEMA(val, ema[i-1], N)
  return ema

We’ll use these functions to calculate the signal line. We’ll buy when the TII crosses above the signal line, and sell/short when the signal line crosses below.

We’ll code this up below:

def TIISignal(data, P=60, N=9, shorts=True):
  df = calcTrendIntensityIndex(data, P=P)
  df['SignalLine'] = calcEMA(df['TII'], N=N)
  df['position'] = np.nan
  df['position'] = np.where(df['TII']>=df['SignalLine'], 1, 
                            df['position'])
  if shorts:
    df['position'] = np.where(df['TII']<df['SignalLine'], -1,
                              df['position'])
  else:
    df['position'] = np.where(df['TII']<df['SignalLine'], 0,
                              df['position'])
    
  df['position'] = df['position'].ffill().fillna(0)
  return calcReturns(df)

This strategy is short and sweet. Just buy on those signal line cross-overs.

P = 60
N = 9
df_sig = TIISignal(data.copy(), P=P, N=N)

fig, ax = plt.subplots(3, figsize=(12, 8), sharex=True)
ax[0].plot(df_sig['Close'], label='Close')
ax[0].plot(df_sig[f'SMA_{P}'], label=f'SMA({P})')
ax[0].set_ylabel('Price ($)')
ax[0].set_title(f'Price and {P}-Day SMA for {ticker}')
ax[0].legend()

ax[1].plot(df_sig['TII'], label='TII')
ax[1].plot(df_sig['SignalLine'], label='Signal Line')
ax[1].set_ylabel('TII')
ax[1].set_title('Trend Intensity Index and Signal Line')
ax[1].legend(bbox_to_anchor=[1, 0.65])

ax[2].plot(df_sig['strat_cum_returns']*100, label='Signal Line')
ax[2].plot(df_sig['cum_returns']*100, label='Buy and Hold')
ax[2].set_title('Cumulative Returns')
ax[2].set_ylabel('Returns (%)')
ax[2].set_xlabel('Date')
ax[2].legend()
plt.tight_layout()
plt.show()

df_stats = pd.concat([df_stats,
                      pd.DataFrame(getStratStats(
                        df_sig['strat_log_returns']),
                        index=['Signal'])])
df_stats

This model missed out on the big gold bull market from 2008-2011, enduring a long and protracted drawdown. It did have some great returns as the gold price fell into 2016. After running flat for a few years, closely following the GDX, it started to pick up with the underlying and managed to get on the right side of the COVID crash – but then gave it all back as the GDX rebounded quickly.

Lots of volatility here and not much upside leads to some poor risk adjusted metrics. You can always adjust the standard P and N values used to see if you can find some better parameters for this one.

Trading with Signal Momentum

I decided to throw in a bonus model now that we have the signal line. We can use that as a classic momentum indicator buying when the TII signal crosses the centerline. In this case, I didn’t bother with trying to hit the tops/bottoms when momentum starts to wane like we did above, just buy and sell based on the centerline, which drastically simplifies the code.

def TIISignalMomentum(data, P=60, N=9, centerline=50,
                      shorts=True):
  df = calcTrendIntensityIndex(data, P=P)
  df['SignalLine'] = calcEMA(df['TII'], N=N)
  df['position'] = np.nan
  df['position'] = np.where(df['SignalLine']>centerline, 1, 
                            df['position'])
  if shorts:
    df['position'] = np.where(df['SignalLine']<centerline, -1,
                              df['position'])
  else:
    df['position'] = np.where(df['SignalLine']<centerline, 0,
                              df['position'])
    
  df['position'] = df['position'].ffill().fillna(0)
  return calcReturns(df)

And to test it on our data:

P = 60
N = 9
df_sigM = TIISignalMomentum(data.copy(), P=P, N=N)

fig, ax = plt.subplots(3, figsize=(12, 8), sharex=True)
ax[0].plot(df_sigM['Close'], label='Close')
ax[0].plot(df_sigM[f'SMA_{P}'], label=f'SMA({P})')
ax[0].set_ylabel('Price ($)')
ax[0].set_title(f'Price and {P}-Day SMA for {ticker}')
ax[0].legend()

ax[1].plot(df_sigM['TII'], label='TII')
ax[1].plot(df_sigM['SignalLine'], label='Signal Line')
ax[1].axhline(50, label='Centerline', c='k', linestyle=':')
ax[1].set_ylabel('TII')
ax[1].set_title('Trend Intensity Index and Signal Line')
ax[1].legend(bbox_to_anchor=[1, 0.65])

ax[2].plot(df_sigM['strat_cum_returns']*100, label='Signal Line')
ax[2].plot(df_sigM['cum_returns']*100, label='Buy and Hold')
ax[2].set_title('Cumulative Returns')
ax[2].set_ylabel('Returns (%)')
ax[2].set_xlabel('Date')
ax[2].legend()
plt.tight_layout()
plt.show()

df_stats = pd.concat([df_stats,
                      pd.DataFrame(getStratStats(
                        df_sigM['strat_log_returns']),
                        index=['Signal Momentum'])])
df_stats

This wound up being the worst performer of them all! This strategy lost over 70% of its initial capital and spent 17 years in a drawdown. I don’t know anyone who would have the patience to stick with a strategy to endure those kinds of losses.

Of course, you’ve got the standard knobs to play with to try to pull some better returns out of this one. Plus you could try adding the complexity from the previous momentum strategy to see how that could boost returns.

Improving Your Trading

The TII is designed to measure the strength of a trend. We walked through four different trading strategies, but none really shined. That could be due to a variety of factors – poor selection of the underlying security, bad parameterization, it ought to be combined with other indicators, or maybe this indicator is just terrible.

The only way to know for sure is by testing it.

Proper testing and trading is difficult. The code here should not be relied upon for running your account – it makes too many simplifying assumptions for a real system. It’s just here to teach you about various indicators, running a backtest, and give you a feel for some of the key metrics used to evaluate a strategy.

If you want a fuller backtest experience, check us out where we provide professional backtests and data in a no-code framework so you can test ideas and find valuable trading ideas quickly.

Should you Trade with the Kelly Criterion?

The Kelly Criterion gives an optimal result for betting based on the probability of winning a bet and how much you receive for winning. If you check out Wikipedia or Investopedia, you’ll see formulas like this:

f^{*} = p - \frac{1-p}{b-1}

which gives you the optimal amount to bet (f^*) given the probability of winning (p) and the payout you’re given for the bet (b). For example, if you have a bet that gives you a 52% chance of winning and you have 1:1 odds (bet 1 to win 1 and your money back, so you get a total payout of 2), then you should bet 4% of your capital on each game (0.52 – 0.48/(2-1) = 0.04).

This is fine for a binary, win-lose outcome. The trouble is, investing in stocks don’t follow this simple model. If you make a winning trade, you could get a 10% return, 8% return, 6.23%, 214%, or any other value.
So how do we change our binary formula to a continuous model?

Continuous Kelly

Ed Thorpe and Claude Shannon (yes, the Claude Shannon for us nerds out there) used the Kelly Criterion to manage their black jack bankroll and clean up in Vegas in the 60s. Seeing the applicability of this method, Thorpe extended it to Wall Street using it to manage his investments while running his own hedge fund. Developing a continuous Kelly model isn’t actually very straightforward, but Thorpe offers this series of steps which you can find in Section 7 here.

Most people probably don’t care about the derivation, so we’ll just jump to this new model which winds up being deceptively simple:

f^{*} = \frac{\mu - r}{\sigma^2}

Here, \mu are the mean returns, r is the risk-free rate, and \sigma is the standard deviation of returns. All of this looks very much like the Sharpe Ratio, but instead of dividing by the standard deviation, we divide by the variance.

Now that we have our new formula, let’s put an example strategy or two together to see how it works in theory.

Long-Only Portfolio with the Kelly Criterion

We’re going to start simple to put this into practice with a long-only portfolio that will rebalance based on the Kelly Criterion. When trading with Kelly position sizing, there are a few things to keep in mind.

First, our Kelly factor (f^*) can go over 1, which implies use of leverage. This may or may not be feasible depending on the trader and leverage may incur additional borrowing costs. Clearly, it increases your risk and volatility as well so not everyone will want to trade with leverage.

Second, we need to estimate the mean and standard deviation of our security. This can be affected by our look-back period. Do we use one year? One month? Or some other time period? One year seems to be standard, so we’ll start with that. Also, the time period you choose may be related to the speed of your trading strategy. If you’re trading minute-by-minute, perhaps a one-month look-back period makes more sense.

Finally, we need to determine how frequently we’ll rebalance our portfolio according to our updated Kelly factor. This can create a lot of transaction costs if it occurs too frequently. There could also be tax implications associated with selling positions quickly, which are beyond the scope of this model. For our purposes, we’ll rebalance every day so that we’re always holding the estimated, optimal position.

Ok, with that out of the way, let’s get some Python packages.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import yfinance as yf

Calculating the Kelly Criterion for Stocks

Our first function is going to be straightforward. All we’re doing is plugging the mean, standard deviation, and our interest rate into that formula above and returning f-star.

def calcKelly(mean, std, r):
  return (mean - r) / std**2

Our example is going to be a simple, vectorized backtest to quickly give you the flavor of using the Kelly formula for position sizing. We’ll write one more helper function that will take our returns as a Pandas Series and return the f-star for each time step.

def getKellyFactor(returns: pd.Series, r=0.01, 
  max_leverage=None, periods=252, rolling=True):
  '''
  Calculates the Kelly Factor for each time step based
  on the parameters provided.
  '''
  if rolling:
    std = returns.rolling(periods).std()
    mean = returns.rolling(periods).mean()
  else:
    std = returns.expanding(periods).std()
    mean = returns.expanding(periods).mean()
  r_daily = np.log((1 + r) ** (1 / 252))
  kelly_factor = calcKelly(mean, std, r_daily)
  # No shorts
  kelly_factor = np.where(kelly_factor<0, 0, kelly_factor)
  if max_leverage is not None:
    kelly_factor = np.where(kelly_factor>max_leverage,
      max_leverage, kelly_factor)
    
  return kelly_factor

Let’s use some real data to demonstrate this. I chose the SPY, which is the S&P 500 ETF which was introduced back in 1993. We can get it’s history through 2020 from YFinance with the code below:

ticker = 'SPY'
yfObj = yf.Ticker(ticker)
data = yfObj.history(start='1993-01-01', end='2020-12-31')
# Drop unused columns
data.drop(['Open', 'High', 'Low', 'Volume', 'Dividends', 
  'Stock Splits'], axis=1, inplace=True)
data.head()

Now we’re ready to put this into action in a strategy. A key assumption stated by Thorpe in the Kelly formula above is the lack of short selling. With that in mind, we’re going to start with the simplest strategy I can think of, a buy-and-hold strategy that will rebalance the portfolio between cash and equities according the Kelly factor.

For this, we’re going to assign a percentage of our total portfolio to the SPY based on the Kelly factor. If leverage increases above 1, then we will have a negative cash value to reflect our borrowing. I didn’t get very granular, so I just used the same interest rate from 1993-2020 and assumed you can borrow and lend at that same rate – which is really a ridiculous assumption because nobody is going to let you, dear retail investor, pay 1% to lever up your stock portfolio (commodities, FOREX, crypto, and other markets typically have more leverage available). You can get daily 10-year rates or whatever treasury instrument suits you as a baseline from FRED, the US Treasury site, or your favorite data provider to use real data in your Kelly calculation.

Let’s get to the function.

def LongOnlyKellyStrategy(data, r=0.02, max_leverage=None, periods=252, 
  rolling=True):
  data['returns'] = data['Close'] / data['Close'].shift(1)
  data['log_returns'] = np.log(data['returns'])
  data['kelly_factor'] = getKellyFactor(data['log_returns'], 
    r, max_leverage, periods, rolling)
  cash = np.zeros(data.shape[0])
  equity = np.zeros(data.shape[0])
  portfolio = cash.copy()
  portfolio[0] = 1
  cash[0] = 1
  for i, _row in enumerate(data.iterrows()):
    row = _row[1]
    if np.isnan(row['kelly_factor']):
      portfolio[i] += portfolio[i-1]
      cash[i] += cash[i-1]
      continue

    portfolio[i] += cash[i-1] * (1 + r)**(1/252) + equity[i-1] * row['returns']
    equity[i] += portfolio[i] * row['kelly_factor']
    cash[i] += portfolio[i] * (1 - row['kelly_factor'])

  data['cash'] = cash
  data['equity'] = equity
  data['portfolio'] = portfolio
  data['strat_returns'] = data['portfolio'] / data['portfolio'].shift(1)
  data['strat_log_returns'] = np.log(data['strat_returns'])
  data['strat_cum_returns'] = data['strat_log_returns'].cumsum()
  data['cum_returns'] = data['log_returns'].cumsum()
  return data

We’re going to run a max-leverage strategy and compare it to a baseline, buy-and-hold.

kelly = LongOnlyKellyStrategy(data.copy())

fig, ax = plt.subplots(2, figsize=(12, 8), sharex=True)
ax[0].plot(np.exp(kelly['cum_returns']) * 100, label='Buy and Hold')
ax[0].plot(np.exp(kelly['strat_cum_returns']) * 100, label='Kelly Model')
ax[0].set_ylabel('Returns (%)')
ax[0].set_title('Buy-and-hold and Long-Only Strategy with Kelly Sizing')
ax[0].legend()

ax[1].plot(kelly['kelly_factor'])
ax[1].set_ylabel('Leverage')
ax[1].set_xlabel('Date')
ax[1].set_title('Kelly Factor')

plt.tight_layout()
plt.show()

Our model blew up spectacularly! Without any constraints on our leverage, it quickly shot up to 50x leverage – and did great for a bit – but then hit some major losses and plummeted below the baseline. It wound up losing our initial investment fairly quickly.

Important lesson here, while the Kelly Criterion may give you the optimal allocation to trade with, it is only as good as the assumptions that underpin it. Note that this is *not* a predictive tool, we’re looking back in time and assuming the previous year’s volatility is a good guide going forward.

We can do a few different things to de-risk this while gaining the benefits of a better money management strategy. First, we can cap our leverage to something more reasonable, say 3 or 4x (although I’d even keep it lower myself). Let’s see how this impacts our model.

def getStratStats(log_returns: pd.Series,
  risk_free_rate: float = 0.02):
  stats = {}  # Total Returns
  stats['tot_returns'] = np.exp(log_returns.sum()) - 1  
  
  # Mean Annual Returns
  stats['annual_returns'] = np.exp(log_returns.mean() * 252) - 1  
  
  # Annual Volatility
  stats['annual_volatility'] = log_returns.std() * np.sqrt(252)  
  
  # Sortino Ratio
  annualized_downside = log_returns.loc[log_returns<0].std() * \
    np.sqrt(252)
  stats['sortino_ratio'] = (stats['annual_returns'] - \
    risk_free_rate) / annualized_downside  
  
  # Sharpe Ratio
  stats['sharpe_ratio'] = (stats['annual_returns'] - \
    risk_free_rate) / stats['annual_volatility']  
  
  # Max Drawdown
  cum_returns = log_returns.cumsum() - 1
  peak = cum_returns.cummax()
  drawdown = peak - cum_returns
  max_idx = drawdown.argmax()
  stats['max_drawdown'] = 1 - np.exp(cum_returns[max_idx]) / np.exp(peak[max_idx])
  
  # Max Drawdown Duration
  strat_dd = drawdown[drawdown==0]
  strat_dd_diff = strat_dd.index[1:] - strat_dd.index[:-1]
  strat_dd_days = strat_dd_diff.map(lambda x: x.days).values
  strat_dd_days = np.hstack([strat_dd_days,
    (drawdown.index[-1] - strat_dd.index[-1]).days])
  stats['max_drawdown_duration'] = strat_dd_days.max()
  return stats

max_leverage = np.arange(1, 6)

fig, ax = plt.subplots(2, figsize=(15, 10), sharex=True)
data_dict = {}
df_stats = pd.DataFrame()
for l in max_leverage:
  kelly = LongOnlyKellyStrategy(data.copy(), max_leverage=l)
  data_dict[l] = kelly.copy()

  ax[0].plot(np.exp(kelly['strat_cum_returns']) * 100,
             label=f'Max Leverage = {l}')

  ax[1].plot(kelly['kelly_factor'], label=f'Max Leverage = {l}')

  stats = getStratStats(kelly['strat_log_returns'])
  df_stats = pd.concat([df_stats, 
    pd.DataFrame(stats, index=[f'Leverage={l}'])])

ax[0].plot(np.exp(kelly['cum_returns']) * 100, label='Buy and Hold', 
           linestyle=':')
ax[0].set_ylabel('Returns (%)')
ax[0].set_title('Buy-and-hold and Long-Only Strategy with Kelly Sizing')
ax[0].legend()

ax[1].set_ylabel('Leverage')
ax[1].set_xlabel('Date')
ax[1].set_title('Kelly Factor')

plt.tight_layout()
plt.show()

stats = pd.DataFrame(getStratStats(kelly['log_returns']), index=['Buy and Hold'])
df_stats = pd.concat([stats, df_stats])
df_stats

There’s a lot to unpack here depending on the leverage ratio. For clarity, a leverage ratio of 1, means that we aren’t actually using any leverage, we’re maxing out by putting all of our portfolio into the SPY. This more conservative model, winds up with lower total returns than the buy and hold approach, but it has a shorter drawdown with less volatility.

Allowing leverage to max out at 2x increases the total returns over buy-and-hold, but with more volatility and lower risk-adjusted returns. It’s a bit tough to see these models in the plot above, so we’ll zoom in on the these results below.

fig, ax = plt.subplots(2, figsize=(15, 8), sharex=True)

ax[0].plot(np.exp(data_dict[1]['strat_cum_returns'])*100, 
         label='Max Leverage=1')
ax[0].plot(np.exp(data_dict[2]['strat_cum_returns'])*100, 
         label='Max Leverage=2')
ax[0].plot(np.exp(data_dict[1]['cum_returns'])*100, 
         label='Buy and Hold', linestyle=':')
ax[0].set_ylabel('Returns (%)')
ax[0].set_title('Low-Leverage and Baseline Models')
ax[0].legend()

ax[1].plot(data_dict[1]['kelly_factor'], label='Max Leverage = 1')
ax[1].plot(data_dict[2]['kelly_factor'], label='Max Leverage = 2')
ax[1].set_xlabel('Date')
ax[1].set_ylabel('Leverage')
ax[1].set_title('Kelly Factor')

plt.tight_layout()
plt.show()

Here, we can see that the Kelly Criterion tends to get out of the market and go to cash as the volatility increases during large drawdowns. The worst of the crashes in 2000 and 2008 are avoided. The COVID crash in 2020, however, was much more rapid and wound up leading to the big losses – particularly the more leveraged the strategy was – as the strategy stayed in the market during the worst of it and got crushed.

Unsurprisingly, the more highly leveraged models all had more volatility with much larger drawdowns. They gave us higher total returns, but the risk metrics all show they did so with a lot more risk.

Above, I mentioned a few ways we can cut down on our leverage beyond just a hard cap.

Larger Sample Size

We can increase the sample size by increasing the number of periods we use for calculating the mean and standard deviation of our returns. This will lead to a slower reacting model, but may lead to better estimates and thus better results.

max_leverage = 3

periods = 252 * np.arange(1, 5)

fig, ax = plt.subplots(2, figsize=(15, 10), sharex=True)
data_dict = {}
df_stats = pd.DataFrame()
for p in periods:
  p = int(p)
  kelly = LongOnlyKellyStrategy(data.copy(), periods=p,
      max_leverage=max_leverage)
  data_dict[p] = kelly.copy()

  ax[0].plot(np.exp(kelly['strat_cum_returns']) * 100,
             label=f'Days = {p}')

  ax[1].plot(kelly['kelly_factor'], label=f'Days = {p}', linewidth=0.5)

  stats = getStratStats(kelly['strat_log_returns'])
  df_stats = pd.concat([df_stats, 
    pd.DataFrame(stats, index=[f'Days={p}'])])

ax[0].plot(np.exp(kelly['cum_returns']) * 100, label='Buy and Hold', 
           linestyle=':')
ax[0].set_ylabel('Returns (%)')
ax[0].set_title(
    'Buy-and-hold and Long-Only Strategy with Kelly Sizing ' +
    'and Variable Lookback Periods')
ax[0].legend()

ax[1].set_ylabel('Leverage')
ax[1].set_xlabel('Date')
ax[1].set_title('Kelly Factor')

plt.tight_layout()
plt.show()

stats = pd.DataFrame(getStratStats(kelly['log_returns']), index=['Buy and Hold'])
df_stats = pd.concat([stats, df_stats])
df_stats

These longer lookback periods tend to increase the time out of market, which extends the drawdown durations in many cases. But they wind up with higher overall returns, at least up to the 4-year lookback period, One thing that makes these comparisons somewhat unfair is that they require more time and data before they enter the market in the first place. So the SPY gets to compounding immediately, the 1-year model gets going after 252 trading days have passed, and so forth. Regardless, we still see good gains from these longer term systems.

One more example of this before moving on. I included an argument called rolling in our function that defaults to True. It’s almost always better to include more data in our estimates, so if we set this to False, we switch out a rolling horizon for an expanding horizon calculation. This latter approach will be identical to the rolling horizon model at day 252 under the standard settings, but then will diverge because it doesn’t drop old data. So the mean and standard deviation of the SPY will contain the last 500 data points at day 500. This quickly increases our sample size as the model moves forward in time. In addition, the number of periods we use now serves as a minimum number of data points required before we get a result for our Kelly factor.

kelly = LongOnlyKellyStrategy(data.copy(), rolling=False)

fig, ax = plt.subplots(2, figsize=(12, 8), sharex=True)
ax[0].plot(np.exp(kelly['cum_returns']) * 100, label='Buy and Hold')
ax[0].plot(np.exp(kelly['strat_cum_returns']) * 100, label='Kelly Model')
ax[0].set_ylabel('Returns (%)')
ax[0].set_title(
    'Buy-and-hold and Long-Only Strategy with Kelly Sizing ' +
    'and Expanding Horizon')
ax[0].legend()

ax[1].plot(kelly['kelly_factor'])
ax[1].set_ylabel('Leverage')
ax[1].set_xlabel('Date')
ax[1].set_title('Expanding Kelly Factor')

plt.tight_layout()
plt.show()

In this case, we see that the Kelly factor largely decreases over time as the sample grows. It does begin to increase after 2008, however, it performed very poorly after the 2000 crash – which it was highly levered heading into – and was never able to recover.

Kelly Sizing and a Trading Strategy

So far we have just looked at applying the Kelly Criterion to a single asset to manage our cash-equity balance. What if we want to use it with a trading strategy on a single asset? Could we improve our risk adjusted returns in this scenario.
The answer is “yes,” but we do need to be careful in how we apply it. If we run the risk of blowing up thanks to leverage while trading the S&P 500, this can be even more of a risk with a given trading strategy.

We’re going to keep things as simple as possible here and run a moving average cross-over strategy. Again, we’ll keep all the assumptions about liquidity, leverage, and no short selling that we laid out above.

There are a couple of changes we’ll need to make to run this. First, we will need to figure out our position, which is just going to be when the short-term SMA is above the long-term SMA.

Second, instead of using the mean and standard deviation of the underlying asset, we’re going to rely on the stats from our strategy. This will use the stats from using our strategy without position sizing.

The code is given below and is similar to what you saw in our long-only strategy.

# Kelly money management for trading strategy
def KellySMACrossOver(data, SMA1=50, SMA2=200, r=0.01, 
  periods=252, max_leverage=None, rolling=True):
  '''
  Sizes a simple moving average cross-over strategy according
  to the Kelly Criterion.
  '''
  data['returns'] = data['Close'] / data['Close'].shift(1)
  data['log_returns'] = np.log(data['returns'])
  # Calculate positions
  data['SMA1'] = data['Close'].rolling(SMA1).mean()
  data['SMA2'] = data['Close'].rolling(SMA2).mean()
  data['position'] = np.nan
  data['position'] = np.where(data['SMA1']>data['SMA2'], 1, 0)
  data['position'] = data['position'].ffill().fillna(0)

  data['_strat_returns'] = data['position'].shift(1) * \
    data['returns']
  data['_strat_log_returns'] = data['position'].shift(1) * \
    data['log_returns']

  # Calculate Kelly Factor using the strategy's returns
  kf = getKellyFactor(data['_strat_log_returns'], r, 
    max_leverage, periods, rolling)
  data['kelly_factor'] = kf

  cash = np.zeros(data.shape[0])
  equity = np.zeros(data.shape[0])
  portfolio = cash.copy()
  portfolio[0] = 1
  cash[0] = 1
  for i, _row in enumerate(data.iterrows()):
    row = _row[1]
    if np.isnan(kf[i]):
      portfolio[i] += portfolio[i-1]
      cash[i] += cash[i-1]
      continue

    portfolio[i] += cash[i-1] * (1 + r)**(1/252) + equity[i-1] * row['returns']
    equity[i] += portfolio[i] * row['kelly_factor']
    cash[i] += portfolio[i] * (1 - row['kelly_factor'])

  data['cash'] = cash
  data['equity'] = equity
  data['portfolio'] = portfolio
  data['strat_returns'] = data['portfolio'] / data['portfolio'].shift(1)
  data['strat_log_returns'] = np.log(data['strat_returns'])
  data['strat_cum_returns'] = data['strat_log_returns'].cumsum()
  data['cum_returns'] = data['log_returns'].cumsum()
  return data

Let’s see how the model performs vs the baseline with moderate leverage.

kelly_sma = KellySMACrossOver(data.copy(), max_leverage=3)

fig, ax = plt.subplots(2, figsize=(15, 8), sharex=True)
ax[0].plot(np.exp(kelly_sma['cum_returns']) * 100, label='Buy-and-Hold')
ax[0].plot(np.exp(kelly_sma['strat_cum_returns'])* 100, label='SMA-Kelly')
ax[0].plot(np.exp(kelly_sma['_strat_log_returns'].cumsum()) * 100, label='SMA')
ax[0].set_ylabel('Returns (%)')
ax[0].set_title('Moving Average Cross-Over Strategy with Kelly Sizing')
ax[0].legend()

ax[1].plot(kelly_sma['kelly_factor'])
ax[1].set_ylabel('Leverage')
ax[1].set_xlabel('Date')
ax[1].set_title('Kelly Factor')

plt.tight_layout()
plt.show()

sma_stats = pd.DataFrame(getStratStats(kelly_sma['log_returns']), 
                         index=['Buy and Hold'])
sma_stats = pd.concat([sma_stats,
            pd.DataFrame(getStratStats(kelly_sma['strat_log_returns']),
              index=['Kelly SMA Model'])])
sma_stats = pd.concat([sma_stats,
            pd.DataFrame(getStratStats(kelly_sma['_strat_log_returns']),
              index=['SMA Model'])])
sma_stats

The Kelly SMA Model doubles the buy and hold approach in terms of total returns. Like the others, it was in a leveraged long position heading into the COVID crash and got crushed. Moreover, from a risk-adjusted return basis, it performs worse than either a basic SMA model or the buy and hold strategy.

Trading with the Kelly Criterion

Leverage can be a powerful and dangerous tool by amplifying both gains and losses. Each of the strategies we ran without a cap on the leverage ratio blew up at some point, highlighting the dangers of such an approach. Most traders who do use the Kelly Criterion in their position sizing only trade half or quarter Kelly, i.e. with 50% or 25% of the Kelly factor size. This is to control risk and avoid blowing up, which is a fate much worse than underperforming the market for a few years.

It’s important to note that the Kelly Criterion is not predictive. While it may optimize the long-run growth of your returns, we see that it falls down time and time again when unconstrained because it is looking backwards over a small sample size. To get the most out of it, we’d need to use it in a constrained setting on a diversified strategy over many markets. That way, we’d be developing our stats based on the performance of the strategy, giving us a larger sample size and better estimate while also limiting our leverage. We’ll look at strategies like this in future posts.

For now, if you want to keep up with what we’re doing at Raposa, sign up with your email below!

How to Trade with the MACD

The Moving Average Convergence-Divergence (MACD, sometimes pronounced “Mac-Dee”) is a very popular indicator that is frequently used in momentum and trend following systems. In short, the MACD is a trailing indicator that gives an indication of the general price trend of a security.

TL;DR

We walk through the reasoning and math behind the MACD along with Python code so you can learn to apply it yourself.

How the MACD Works

Despite the intimidating name, the MACD is relatively straightforward to understand and calculation. It takes some ideas from other indicators, namely EMA and moving average cross-over strategies, and combines them into a single, easy to use value.

In its most basic form, we have two EMA signals, a fast one and a slow one (12 and 26-days are popularly chosen). These are calculated every day, then we subtract the fast one from the slow one to get the difference. In psuedocode we have:

  1. Calculate fast EMA:
    • EMA_fast[t] = (Price - EMA_fast[t-1]) * 2 / (N_fast + 1) + EMA_fast[t-1] 
  2. Calculate slow EMA:
    • EMA_slow[t] = (Price - EMA_slow[t-1]) * 2 / (N_slow + 1) + EMA_slow[t-1]
  3. Subtract the two to get MACD:
    • MACD[t] = EMA_fast[t] - EMA_slow[t]

Or, if you prefer to write it mathematically:

EMA_t = \big( P_t - EMA_{t-1} \big) \frac{2}{N+1} + EMA_{t-1}
MACD_t = EMA_t^{fast} - EMA_t^{slow}

It’s really that easy.

The MACD is called the Moving Average Convergence-Divergence. What are converging and diverging, are the two EMAs. As the short-term EMA and long-term EMA converge, they get closer to the same value and the indicator moves to 0. Divergence is driven by the short-term EMA moving up or down causing the distance between the two EMAs to move farther apart. There are other ways the terms convergence and divergence are used with this indicator which we’ll get to below.

The most basic way to trade it is buy when MACD > 0, and sell/short when MACD < 0. When it’s positive, we have the faster EMA above the longer EMA, and vice versa when it goes negative. This set up is a basic, exponential moving average crossover strategy.

MACD Example in Python

To code it, we just need a few basic packages.

import numpy as np
import pandas as pd
import yfinance as yf
import matplotlib.pyplot as plt

The details for the EMA calculation are given here, so we’ll just show the final code which we’ll leverage in our model.

def _calcEMA(P, last_ema, N):
    return (P - last_ema) * (2 / (N + 1)) + last_ema

def calcEMA(data: pd.DataFrame, N: int, key: str = 'Close'):
    # Initialize series
    data['SMA_' + str(N)] = data[key].rolling(N).mean()
    ema = np.zeros(len(data)) + np.nan
    for i, _row in enumerate(data.iterrows()):
        row = _row[1]
        if np.isnan(ema[i-1]):
            ema[i] = row['SMA_' + str(N)]
        else:
            ema[i] = _calcEMA(row[key], ema[i-1], N)

    data['EMA_' + str(N)] = ema.copy()
    return data

With the calcEMA function above, we can easily write our MACD function. We just need to call calcEMA twice with the fast and slow parameters, and subtract the values.

def calcMACD(data: pd.DataFrame, N_fast: int, N_slow: int):
    assert N_fast < N_slow, "Fast EMA must be less than slow EMA parameter."
    # Add short term EMA
    data = calcEMA(data, N_fast)
    # Add long term EMA
    data = calcEMA(data, N_slow)
    # Subtract values to get MACD
    data['MACD'] = data[f'EMA_{N_fast}'] - data[f'EMA_{N_slow}']
    return data

We’re ready to test it. I just grabbed a one year time period from a random stock in the S&P 500 to illustrate how the indicator works. We’ll just run it through our function and take a look at the output.

ticker = 'DTE'
start = '2013-01-01'
end = '2014-01-01'
yfObj = yf.Ticker(ticker)
df = yfObj.history(start=start, end=end)
N_fast = 12
N_slow = 26
data = calcMACD(df, N_fast, N_slow)
# Drop extra columns
data.drop(['Open', 'High', 'Low', 'Volume', 'Dividends', 'Stock Splits'], axis=1, inplace=True)
data.iloc[N_slow-5:N_slow+5]

As you can see in the table above, our function provides the different EMA values (and the SMAs they’re initialized on) in addition to the MACD. We can plot these values below to see how they track with one another.

fig, ax = plt.subplots(2, figsize=(12, 8), sharex=True)
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

ax[0].plot(data['Close'], label=f'{ticker}', linestyle=':')
ax[0].plot(data[f'EMA_{N_fast}'], label=f'EMA-{N_fast}')
ax[0].plot(data[f'EMA_{N_slow}'], label=f'EMA-{N_slow}')
ax[0].set_ylabel('Price ($)')
ax[0].set_title(f'Price and EMA Values for {ticker}')
ax[0].legend()

ax[1].plot(data['MACD'], label='MACD')
ax[1].set_ylabel('MACD')
ax[1].set_xlabel('Date')
ax[1].set_title(f'MACD for {ticker}')

plt.show()

We have a few good oscillations around 0, which, in our simplest strategy, would indicate buy/sell signals as the EMAs cross over. You can also see that the initial price rise from January through May has a high MACD (>0.5), and the MACD peak roughly corresponds with the peak in the price of our security. The MACD begins to fall and turn negative, which follows a healthy retracement in the price before it rebounds again.

Also, note that as the MACD becomes more positive, this means that the short-term EMA is increasing more rapidly than the longer, slower EMA we’re comparing it to. This is a sign of stronger upward momentum in the price action. The opposite holds true when the MACD becomes more negative.

There’s another way to generate signals from the MACD, however. This method is based on using a signal line in addition to the MACD as we calculated it.

The MACD with a Signal Line

The signal line is a EMA of the MACD signal we calculated. Frequently this will be a 9-day EMA to go along with a 12 and 26-day EMA like we calculated above.

Writing the psuedocode, we just need to add one more step:

  1. Calculate the Signal Line (SL) as the EMA of the MACD:
    • SL[t] = (MACD[t] - EMA_MACD[t-1]) * 2 / (N_SL + 1) + EMA_MACD[t-1]

Again, we can write this mathematically as:

SL_t = \big( MACD_t - EMA_{t-1}^{MACD} \big) \frac{2}{N^{SL}+1} + EMA_{t-1}^{MACD}

Let’s modify our calcMACD function to allow for signal line calculation too.

def calcMACD(data: pd.DataFrame, N_fast: int, N_slow: int, 
             signal_line: bool = True, N_sl: int = 9):
    assert N_fast < N_slow, "Fast EMA must be less than slow EMA parameter."
    # Add short term EMA
    data = calcEMA(data, N_fast)
    # Add long term EMA
    data = calcEMA(data, N_slow)
    # Subtract values to get MACD
    data['MACD'] = data[f'EMA_{N_fast}'] - data[f'EMA_{N_slow}']
    if signal_line:
        data = calcEMA(data, N_sl, key='MACD')
        # Rename columns
        data.rename(
            columns={f'SMA_{N_sl}': f'SMA_MACD_{N_sl}',
                     f'EMA_{N_sl}': f'SignalLine_{N_sl}'},
           inplace=True)
    return data

Now, we will run this new function and plot it to show how the signal line looks versus the MACD.

N_fast = 12
N_slow = 26
N_sl = 9
signal_line = True
data = calcMACD(df, N_fast, N_slow, signal_line, N_sl)

# Plot MACD and Signal Line
fig, ax = plt.subplots(figsize=(15, 8))
ax.plot(data['MACD'], label='MACD')
ax.plot(data[f'SignalLine_{N_sl}'], label='Signal Line')
ax.set_ylabel('MACD')
ax.set_xlabel('Date')
ax.set_title(f'MACD and Signal Line for {ticker}')
ax.legend()
plt.show()

In the plot above, we can see that the signal line intersects with the MACD on a number of occasions. As you have probably guessed, these intersections provide buy and sell signals as well. You can buy when the MACD crosses above the signal line and sell/short when it goes below.

MACD Momentum

Often times, you’ll see MACD charts with bars as well as lines on the graph like this:

data['MACD_bars'] = data['MACD'] - data[f'SignalLine_{N_sl}']
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

fig, ax = plt.subplots(1, figsize=(15, 8), sharex=True)

# Re-index dates to avoid gaps over weekends
_x_axis = np.arange(data.shape[0])
month = pd.Series(data.index.map(lambda x: x.month))
x = month - month.shift(1)
_x_idx = np.where((x.values!=0) &
                  (~np.isnan(x.values)))[0]

ax.plot(_x_axis, data['MACD'], label='MACD')
ax.plot(_x_axis, data[f'SignalLine_{N_sl}'], label='Signal Line',
       c=colors[2])
ax.bar(_x_axis, data['MACD_bars'], label='MACD Bars', 
      color=colors[4], width=1, edgecolor='w')
ax.set_xticks(_x_axis[np.where(~np.isnan(_x_dates))])
ax.set_xticks(_x_idx)
ax.set_xticklabels(data.index[_x_idx].map(
    lambda x: datetime.strftime(x, '%Y-%m-%d')), 
    fontsize=10)

ax.legend()
plt.show()

The MACD bars, also called MACD histogram, plots the difference between the MACD and the signal line. These bars can be used by chartists (those who read charts to make trades) to determine the strength of the momentum and make buy sell decisions. If the bars are growing in size, then momentum is increasing as the MACD is pulling away from the signal line. Additionally, if the bars begin shrinking, we could be seeing a reversal coming.

We’re focused on automatic and systematic trading, so we don’t want to spend time looking at bars to see if they’re growing or shrinking. Rather, we want the computer to do it for us!

Thankfully, this is only a few small steps away.

Getting Momentum from the MACD Histogram

First, we need to determine how many consecutive days we want before we determine whether or not we want to enter a position. For our illustration here, we’ll use three days. From there, we can get the daily changes by using the diff() method and then use the np.sign() function to determine if it’s positive or negative (note that this function gives 0s a positive sign). We’ll call this column Growth.

After that, we use Pandas’ handy rolling() method to sum our Growth column. This will give us 3 if we have three consecutive up days, -3 if we have three straight down days, or some other value between -2 and 2, making it easy for us to pick out our signals. For lack of a better name, let’s call this column Consecutive_Bars.

Finally, we’ll get our position by looking for all of those +/-3 values. We want to go long on increasing upward momentum and short on decreasing downward momentum, so we just multiply the sign of our Consecutive_Bars column by 1 if the absolute value of Consecutive_Bars equals the number of consecutive days.

The code for all of this is given below.

N_consecutive = 3
data['Growth'] = np.sign(data['MACD_bars'].diff(1))
data['Consecutive_Bars'] = data['Growth'].rolling(N_consecutive).sum()
data['Position'] = data['Consecutive_Bars'].map(
    lambda x: np.sign(x) * 1 if np.abs(x) == N_consecutive else 0)

long_idx = data['Position']==1
short_idx = data['Position']==-1

fig, ax = plt.subplots(2, figsize=(15, 8), sharex=True)

ax[0].plot(_x_axis, data['Close'], label=f'{ticker}', linestyle=':')
ax[0].scatter(_x_axis[long_idx], data.loc[long_idx]['Close'], 
              label='Longs', c=colors[3])
ax[0].scatter(_x_axis[short_idx], data.loc[short_idx]['Close'], 
              label='Shorts', c=colors[1])
ax[0].set_ylabel('Price ($)')
ax[0].set_title(f'{ticker} and Long/Short Positions based on MACD Bars')
ax[0].legend()

ax[1].bar(_x_axis, data['MACD_bars'], label='MACD Bars', 
      color=colors[4], width=1, edgecolor='w')
ax[1].scatter(_x_axis[long_idx], data.loc[long_idx]['MACD_bars'], 
              label='Longs', c=colors[3])
ax[1].scatter(_x_axis[short_idx], data.loc[short_idx]['MACD_bars'], 
              label='Shorts', c=colors[1])
ax[1].plot(_x_axis, data['MACD'], label='MACD')
ax[1].plot(_x_axis, data['SignalLine_9'], label='Signal Line',
           color=colors[2])
ax[1].set_ylabel('MACD')
ax[1].set_title('MACD, Signal Line, MACD Bars, and Long/Short Positions')

ax[1].set_xticks(_x_axis[np.where(~np.isnan(_x_dates))])
ax[1].set_xticks(_x_idx)
ax[1].set_xticklabels(data.index[_x_idx].map(
    lambda x: datetime.strftime(x, '%Y-%m-%d')), 
    fontsize=10)

ax[1].legend()
plt.show()

There’s a lot going on in the plots above!

In the first panel, we simply have the price action for our stock and the long/short positions overlaid so you can get a feel for what the stock may be doing when this indicator chooses to go long or short. It starts off catching the end of an uptrend and getting out near the peak, then waits a bit before entering to go short while the stock flat-lines.

As you can see working through the rest of the chart – like any indicator – it has a few mistakes. This indicator also lags the price action by quite a bit. Not only are you working off of three different EMAs, but then you wait for three consecutive days before putting a position on. This allows it to catch some strong trends, but because we were too quick to exit positions (e.g. wait for a single down day in the MACD bars to get out of a long position) we miss out on a good chunk of the moves, even though we were positionally correct. So, it looks like this could be useful, but we may need a better exit strategy to make the most of it.

In the second panel you see the MACD, signal line, MACD bars, and our longs/shorts. We derive the MACD from the price action and mimics it to some extent. The signal line smooths out some of the MACD moves and looks like a wave that is slightly out of phase. You can see our signal on the bars themselves too.

Some concluding thoughts

The MACD is a very popular indicator because of its flexibility and diversity. Here, we discussed just a handful of ways we can interpret the values to develop trading systems, but there are others out there we’ll address in the future.

While developing systems with the MACD, keep in mind that the indicator itself is unbounded, meaning it has no ceiling or floor. This is unlike other oscillating signals like RSI, which have a maximum and a minimum value. Moreover, the value is dependent on price. While our example ranged from a little over 1 to -1, you could have more expensive securities that produce MACD values in the teens or hundreds. This is important because it becomes very difficult – if not impossible – to compare MACD values directly across different assets.

The 12-26-9 MACD with a signal line is a very popular set up, but it isn’t the only one. You’re free to experiment, but note that the farther the short term and long term EMAs move from one another, the larger MACD values you’re going to get. Not that this is necessarily a problem, but it is something to keep an eye on; you may need to adjust some other parameters in your system to compensate.

There’s a lot to take in with this indicator and a lot of moving parts to keep track of. At Raposa, we make this easy for you. You can pick your stocks, set your parameters, and let us handle all of the details for you. We’ll provide you the stats on your backtests and give live signals when you’re ready to trade it in the real world. Sign up below to learn more.

How to Build Your First Momentum Trading Strategy in Python

Have you ever noticed the tendency for stocks that have risen to continue rising in the future? Likewise, a few down days seem to beget more losses. This is known as momentum and strategies that rely on these patterns are momentum-based strategies.

TL;DR

We develop a basic momentum strategy and test it on GameStop to yield strong returns.

Trading Momentum

Traders and investors have long known about the effects of momentum and have found that these effects appear across a wide variety of markets and time frames. Running these strategies on a single instrument is also known as trend following or time series momentum. We’ll stick with the latter name, and abbreviate it TSM from here on out.

While there are a whole host of ways to run this strategy by combining it with a series of indicators, risk management overlays (a must for a real trading system!), and diversification, we’ll start simple by showing how it works on a single instrument.

To get started, we can turn to Python and some standard libraries.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import yfinance as yf

From here, we can build our basic strategy function that we’ll call TSMStrategy. This will take the log returns of our time series, the period we’re interested in, and a boolean variable that determines whether we will allow short positions or not. The function will then simulate the strategy for us and return the cumulative sum of the performance for further analysis.

def TSMStrategy(returns, period=1, shorts=False):
    if shorts:
        position = returns.rolling(period).mean().map(
            lambda x: -1 if x <= 0 else 1)
    else:
        position = returns.rolling(period).mean().map(
            lambda x: 0 if x <= 0 else 1)
    performance = position.shift(1) * returns
    return performance

It’s a simple function for a simple (but powerful) strategy. The next step requires us to provide it with some data and see how it performs. I’m going to choose a stock that his been in the news a lot lately: GameStop (GME) to see how our model performs.

ticker = 'GME'
yfObj = yf.Ticker(ticker)
data = yfObj.history(start='2000-01-01', end='2020-12-31')

Using the yfinance package, we can get the total public history of GME. It began trading in 2002, but setting the start date to 2000 will allow us to pick up the stock from the beginning without any errors.

To pass this to our strategy, we need to calculate the log returns and provide that to our function.

returns = np.log(data['Close'] / data['Close'].shift(1)).dropna()

The simplest TSM we can implement would require us to purchase the stock if it was up yesterday, and sell if it was down (if we’re holding it, otherwise we just wait). Let’s give it a shot.

performance = TSMStrategy(returns, period=1, shorts=False).dropna()
years = (performance.index.max() - performance.index.min()).days / 365
perf_cum = np.exp(performance.cumsum())
tot = perf_cum[-1] - 1
ann = perf_cum[-1] ** (1 / years) - 1
vol = performance.std() * np.sqrt(252)
rfr = 0.02
sharpe = (ann - rfr) / vol
print(f"1-day TSM Strategy yields:" +
      f"\n\t{tot*100:.2f}% total returns" + 
      f"\n\t{ann*100:.2f}% annual returns" +
      f"\n\t{sharpe:.2f} Sharpe Ratio")
gme_ret = np.exp(returns.cumsum())
b_tot = gme_ret[-1] - 1
b_ann = gme_ret[-1] ** (1 / years) - 1
b_vol = returns.std() * np.sqrt(252)
b_sharpe = (b_ann - rfr) / b_vol
print(f"Baseline Buy-and-Hold Strategy yields:" + 
      f"\n\t{b_tot*100:.2f}% total returns" + 
      f"\n\t{b_ann*100:.2f}% annual returns" +
      f"\n\t{b_sharpe:.2f} Sharpe Ratio")
1-day TSM Strategy yields:
	225.03% total returns
	6.44% annual returns
	0.12 Sharpe Ratio
Baseline Buy-and-Hold Strategy yields:
	184.63% total returns
	5.70% annual returns
	0.07 Sharpe Ratio

The 1-Day TSM strategy beats the buy-and-hold with a reasonable annual return (ignoring transaction costs, which may be high given such a short term strategy). The 1-day lookback is likely fraught with a lot of false trends, so we can run a variety of different time periods to see how they stack up. We’ll just run our model in a loop with 3, 5, 15, 30, and 90-day time periods.

import matplotlib.gridspec as gridspec
periods = [3, 5, 15, 30, 90]
fig = plt.figure(figsize=(12, 10))
gs = fig.add_gridspec(4, 4)
ax0 = fig.add_subplot(gs[:2, :4])
ax1 = fig.add_subplot(gs[2:, :2])
ax2 = fig.add_subplot(gs[2:, 2:])
ax0.plot((np.exp(returns.cumsum()) - 1) * 100, label=ticker, linestyle='-')
perf_dict = {'tot_ret': {'buy_and_hold': (np.exp(returns.sum()) - 1)}}
perf_dict['ann_ret'] = {'buy_and_hold': b_ann}
perf_dict['sharpe'] = {'buy_and_hold': b_sharpe}
for p in periods:
    log_perf = TSMStrategy(returns, period=p, shorts=False)
    perf = np.exp(log_perf.cumsum())
    perf_dict['tot_ret'][p] = (perf[-1] - 1)
    ann = (perf[-1] ** (1/years) - 1)
    perf_dict['ann_ret'][p] = ann
    vol = log_perf.std() * np.sqrt(252)
    perf_dict['sharpe'][p] = (ann - rfr) / vol
    ax0.plot((perf - 1) * 100, label=f'{p}-Day Mean')
    
ax0.set_ylabel('Returns (%)')
ax0.set_xlabel('Date')
ax0.set_title('Cumulative Returns')
ax0.grid()
ax0.legend()
_ = [ax1.bar(i, v * 100) for i, v in enumerate(perf_dict['ann_ret'].values())]
ax1.set_xticks([i for i, k in enumerate(perf_dict['ann_ret'])])
ax1.set_xticklabels([f'{k}-Day Mean' 
    if type(k) is int else ticker for 
    k in perf_dict['ann_ret'].keys()],
    rotation=45)
ax1.grid()
ax1.set_ylabel('Returns (%)')
ax1.set_xlabel('Strategy')
ax1.set_title('Annual Returns')
_ = [ax2.bar(i, v) for i, v in enumerate(perf_dict['sharpe'].values())]
ax2.set_xticks([i for i, k in enumerate(perf_dict['sharpe'])])
ax2.set_xticklabels([f'{k}-Day Mean' 
    if type(k) is int else ticker for 
    k in perf_dict['sharpe'].keys()],
    rotation=45)
ax2.grid()
ax2.set_ylabel('Sharpe Ratio')
ax2.set_xlabel('Strategy')
ax2.set_title('Sharpe Ratio')
plt.tight_layout()
plt.show()

Looking at the above, the 15-day momentum indicator gives us the best absolute and risk-adjusted returns. However, there is a lot of spread in these results. This indicates that we don’t have a very robust strategy (which should be no surprise). We could build on this basic strategy by incorporating other indicators, such as moving averages or exponentially weighted moving averages to indicate momentum. We could manage our risk better by incorporating stop losses or trailing stops to get out of trades closer to the top rather then when we have a down or flat 15-days. We could also allocate capital across multiple securities thereby benefiting from diversification and jumping on trends wherever they appear.

Some of these improvements require a more sophisticated, event-driven backtest system to properly simulate. This is something we’ve been working on for some time at Raposa Technologies and we’re making available to the public. To stay up to date with the latest developments and features as we roll out a professional quality, no-code backtesting system, subscribe below and be among the first to get access and updates to our newest offerings.