Overall Statistics Total Trades2149Average Win0.33%Average Loss-0.32%Compounding Annual Return2.827%Drawdown31.000%Expectancy0.098Net Profit38.460%Sharpe Ratio0.321Loss Rate46%Win Rate54%Profit-Loss Ratio1.02Alpha0.032Beta0Annual Standard Deviation0.1Annual Variance0.01Information Ratio-0.295Tracking Error0.217Treynor Ratio-98.982Total Fees\$0.00
```from QuantConnect.Python import PythonQuandl
from datetime import timedelta
import numpy as np
import pandas as pd

class ImprovedCommodityMomentumTrading(QCAlgorithm):
'''
Demystifying Time-Series Momentum Strategies: Volatility Estimators, Trading Rules and Pairwise Correlations

This paper proposed 3 modifications to the basic time-series momentum strategies in order to reduce portfolio turnover and improve portfolio performance.

1. Volatility Estimator: Yang and Zhang (2000) range-based estimator, which replaces the traditional estimator (standard deviation of past daily returns)
2. Trading Rules: Trading positions takes a continuum of values between -1 and +1 to reflect the statistical strength of price trend, which replaces the traditional trading rules (binary +1 or -1 based on the sign of historical mean return)
3. Pairwise Correlations: Incorporate signed pairwise correlations in the weighing scheme of portfolio construction

Reference:
 Baltas, Nick and Kosowski, Robert, "Demystifying Time-Series Momentum Strategies: Volatility Estimators, Trading Rules and Pairwise Correlations", May 8, 2017.
URL: https://pdfs.semanticscholar.org/a2e9/df201d4b4774fda84a961cc804f2450988c5.pdf
 Yang, Dennis, and Qiang Zhang, "Drift‐Independent Volatility Estimation Based on High, Low, Open, and Close Prices", The Journal of Business, vol. 73, no. 3, 2000, pp. 477–492.
URL: www.jstor.org/stable/10.1086/209650.'''

def Initialize(self):

self.SetStartDate(2008,1, 1)
self.SetEndDate(2019, 9, 1)
self.SetCash(25000)

tickers = ["CHRIS/CME_S1",  # Soybean Futures, Continuous Contract #1
"CHRIS/CME_W1",  # Wheat Futures, Continuous Contract #1
"CHRIS/CME_SM1", # Soybean Meal Futures, Continuous Contract #1
"CHRIS/CME_BO1", # Soybean Oil Futures, Continuous Contract #1
"CHRIS/CME_C1",  # Corn Futures, Continuous Contract #1
"CHRIS/CME_O1",  # Oats Futures, Continuous Contract #1
"CHRIS/CME_LC1", # Live Cattle Futures, Continuous Contract #1
"CHRIS/CME_FC1", # Feeder Cattle Futures, Continuous Contract #1
"CHRIS/CME_LN1", # Lean Hog Futures, Continuous Contract #1
"CHRIS/CME_GC1", # Gold Futures, Continuous Contract #1
"CHRIS/CME_SI1", # Silver Futures, Continuous Contract #1
"CHRIS/CME_PL1", # Platinum Futures, Continuous Contract #1
"CHRIS/ICE_B1",  # Brent Crude Futures, Continuous Contract
"CHRIS/ICE_O1",  # Heating Oil Futures, Continuous Contract #1
"CHRIS/ICE_M1",  # UK Natural Gas Futures, Continuous Contract #1
"CHRIS/ICE_CT1", # Cotton No. 2 Futures, Continuous Contract
"CHRIS/ICE_OJ1", # Orange Juice Futures, Continuous Contract
"CHRIS/ICE_KC1", # Coffee C Futures, Continuous Contract
"CHRIS/ICE_CC1", # Cocoa Futures, Continuous Contract
"CHRIS/ICE_G1",  # Gas Oil Futures, Continuous Contract
"CHRIS/ICE_RS1"] # Canola Futures, Continuous Contract

for ticker in tickers:
data = self.AddData(QuandlFutures, ticker, Resolution.Daily)
data.SetLeverage(3) # Leverage was set to 3 for each of the futures contract

self.OneYear = 365      # time period for trading rule calculation
self.OneMonth = 30      # time period for YZ volatility estimator
self.ThreeMonths = 90   # time period for pairwise correlation calculation

# Last trading date tracker to achieve rebalancing the portfolio every month
self.nextRebalance = self.Time

# Set portfolio target level of volatility, set to 12%
self.portfolio_target_sigma = 0.12

def OnData(self, data):
'''
Monthly rebalance at the beginning of each month.
Portfolio weights for each constituents are calculated based on Baltas and Kosowski weights.
'''

# skip if less than 30 days passed since the last trading date
if self.Time < self.nextRebalance:
return

'''Monthly Rebalance Execution'''
# dataframe that contains the historical data for all securities
history = self.History(self.Securities.Keys, self.OneYear, Resolution.Daily)
history.replace(0, np.nan, inplace = True)

# Get the security symbols are are in the history dataframe
available_symbols = list(set(history.index.get_level_values(level = 0)))

# Liquidate symbols that are not in the history dataframe anymore
for security in self.Securities.Keys:
if security.Value not in available_symbols:
self.Liquidate(security, 'Not found in history request')

# Get the trade signals and YZ volatility for all securities
trade_signals = self.GetTradingSignal(history)
volatility = self.GetYZVolatility(history, available_symbols)

# Get the correlation factor
CF_rho_bar = self.GetCorrelationFactor(history, trade_signals, available_symbols)

#Rebalance the portfolio according to Baltas and Kosowski suggested weights
N_assets = len(available_symbols)
for symbol, signal, vol in zip(available_symbols, trade_signals, volatility):
# Baltas and Kosowski weights (Equation 19 in )
weight = (signal*self.portfolio_target_sigma*CF_rho_bar)/(N_assets*vol)
self.SetHoldings(symbol, weight)

# Set next rebalance time
self.nextRebalance = Expiry.EndOfMonth(self.Time)

def GetCorrelationFactor(self, history, trade_signals, available_symbols):
'''
Calculate the Correlation Factor, which is a function of the average pairwise correlation of all portfolio contituents
- the calculation is based on past three month pairwise correlation
- Notations:
rho_bar - average pairwise correlation of all portfolio constituents
CF_rho_bar - the correlation factor as a function of rho_bar'''

# Get the past three month simple daily returns for all securities
settle = history.settle.unstack(level = 0)
past_three_month_returns = settle.pct_change().loc[settle.index[-1]-timedelta(self.ThreeMonths):]

# Get number of assets
N_assets = len(available_symbols)

# Get the pairwise signed correlation matrix for all assets
correlation_matrix = past_three_month_returns.corr()

# Calculate rho_bar
summation = 0
for i in range(N_assets-1):
for temp in range(N_assets - 1 - i):
j = i + temp + 1
x_i = trade_signals[i]
x_j = trade_signals[j]
rho_i_j = correlation_matrix.iloc[i,j]
summation += x_i * x_j * rho_i_j

# Equation 14 in 
rho_bar = (2 * summation) / (N_assets * (N_assets - 1))

# Calculate the correlation factor (CF_rho_bar)
# Equation 18 in 
return np.sqrt(N_assets / (1 + (N_assets - 1) * rho_bar))

def GetTradingSignal(self, history):
'''
TREND Trading Signal
- Uses the t-statistics of historical daily log-returns to reflect the strength of price movement trend
- TREND Signal Conditions:
t-stat > 1 => TREND Signal = 1
t-stat < 1 => TREND Signal = -1
-1 < t-stat < 1 => TREND Signal = t-stat
'''
settle = history.settle.unstack(level = 0)

# daily futures log-returns based on close-to-close
log_returns = np.log(settle/settle.shift(1)).dropna()

# Calculate the t-statistics as
# (mean-0)/(stdev/sqrt(n)), where n is sample size
mean = np.mean(log_returns)
std = np.std(log_returns)
n = len(log_returns)
t_stat = mean/(std/np.sqrt(n))

# cap holding at 1 and -1
return np.clip(t_stat, a_max=1, a_min=-1)

def GetYZVolatility(self, history, available_symbols):
'''
Yang and Zhang 'Drift-Independent Volatility Estimation'

Formula: sigma_YZ^2 = sigma_OJ^2 + self.k * sigma_SD^2 + (1-self.k)*sigma_RS^2 (Equation 20 in )
where,  sigma_OJ - (Overnight Jump Volitility estimator)
sigma_SD - (Standard Volitility estimator)
sigma_RS - (Rogers and Satchell Range Volatility estimator)'''
YZ_volatility = []

time_index = history.loc[available_symbols].index
today = time_index[-1]

#Calculate YZ volatility for each security and append to list
for ticker in available_symbols:
past_month_ohlc = history.loc[ticker].loc[today-timedelta(self.OneMonth):today]
open, high, low, close = past_month_ohlc.open, past_month_ohlc.high, past_month_ohlc.low, past_month_ohlc.settle
estimation_period = past_month_ohlc.shape

# Calculate constant parameter k for Yang and Zhang volatility estimator
# using the formula found in Yang and Zhang (2000)
k = 0.34 / (1.34 + (estimation_period + 1) / (estimation_period - 1))

# sigma_OJ (overnight jump => stdev of close-to-open log returns)
open_to_close_log_returns = np.log(open/close.shift(1))
open_to_close_log_returns = open_to_close_log_returns[np.isfinite(open_to_close_log_returns)]
sigma_OJ = np.std(open_to_close_log_returns)

# sigma_SD (standard deviation of close-to-close log returns)
close_to_close_log_returns = np.log(close/close.shift(1))
close_to_close_log_returns = close_to_close_log_returns[np.isfinite(close_to_close_log_returns)]
sigma_SD = np.std(close_to_close_log_returns)

# sigma_RS (Rogers and Satchell (1991))
h = np.log(high/open)
l = np.log(low/open)
c = np.log(close/open)
sigma_RS_daily = (h * (h - c) + l * (l - c))**0.5
sigma_RS_daily = sigma_RS_daily[np.isfinite(sigma_RS_daily)]
sigma_RS = np.mean(sigma_RS_daily)

# daily Yang and Zhang volatility
sigma_YZ = np.sqrt(sigma_OJ**2 + k * sigma_SD**2 + (1 - k) * sigma_RS**2)

# append annualized volatility to the list
YZ_volatility.append(sigma_YZ*np.sqrt(252))

return YZ_volatility

class QuandlFutures(PythonQuandl):
def __init__(self):
self.ValueColumnName = "Settle"```