Predicting the daily returns and directions of returns of SPY using technical indicators and fundamental indicators.

#Import the libraries
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
%matplotlib inline
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import AdaBoostRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.covariance import EmpiricalCovariance
from sklearn.preprocessing import StandardScaler #we won't be using this
from sklearn.decomposition import PCA #we won't be using this
from sklearn.linear_model import LogisticRegression 
from sklearn.naive_bayes import GaussianNB
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import AdaBoostClassifier
import xgboost
import sklearn.metrics as metrics
from sklearn.metrics import confusion_matrix
import seaborn as sns
sns.set_style("white")
np.random.seed(7)
from astropy.table import Table, Column

%matplotlib inline

data = pd.read_csv('data.csv', na_values = ".")
data.rename(columns={'Close': 'SPY'}, inplace=True) #renaming the variable close as SPY
data.head(2)

np.array(list(data))  #the csv file we loaded has 28 data series, excluding the index column. The rest of the technical indicators will be calculated along the way. 

#Now that we have read the data, let's see if we have any missing values in any of the columns.
#We will use seaborn's heatmap to do so. 
sns.heatmap(data.isnull(), yticklabels = False, cbar = False, cmap = 'viridis')
#As we can see, there are not missing data. If we had any, there will be yellow marks for the missing rows. 

#Let's make sure that the data types are correct. The period should be a datetime datatype while the rest should be numeric. 
data.ix[:, 1:28] = data.ix[:, 1:28].apply(pd.to_numeric)
data["Period"] = pd.to_datetime(data["Period"])

#Look at the sample size and the first and the last date of the time series. 
print('Sample Size: {}'.format(len(data)))
print(data['Period'].min())
print(data['Period'].max())

data.reset_index(inplace=True)
data_lagged=data.copy()
data_lagged['index'] = data_lagged['index'].apply(lambda x: x+1)
data_merged = data.merge(data_lagged, how="left", on="index", suffixes=("","_prior"))

#Let's now drop columns we won't be needing anymore.
data_merged.drop(['Period_prior', 'Open', 'High', 'Low','Volume', 'VIX', 'MSCI_EM','DP', 'PE', 'BM', 'CAPE',
                  'PCA.price', 'BY', 'DEF', 'TERM', 'CAY','SIM', 'VRP', 'IC', 'BDI', 'NOS', 'CPI', 'PCR',
                  'MA', 'PCA.tech','OIL', 'SI'], axis = 1, inplace = True)
data_merged.dropna(inplace = True)
#This might be a little confusing. The idea is that we should be using whatever we have at a given point 
#in time to predict tomorrow's price. 
#Here, we will be using prior (data at time t) data to predict the return at time t+1. 

#Daily return
rets = data_merged[['SPY', 'SPY_prior', 'VIX_prior','MSCI_EM_prior']]
for i in rets:
    data_merged[i + '_ret'] = np.log(data_merged[i] / data_merged[i].shift(1))

#Trend Indicator
data_merged['SPY_50'] = pd.rolling_mean(data_merged['SPY_prior'], window=50)#this is the 50-day moving average
data_merged['SPY_trend'] = (data_merged['SPY_prior'] - data_merged['SPY_50'])/ data_merged['SPY_50'] 

#Relative Strength Index (RSI)
def relative_strength(prices, n=14):
    deltas = np.diff(prices)
    seed = deltas[:n+1]
    up = seed[seed >= 0].sum()/n
    down = -seed[seed < 0].sum()/n
    rs = up/down
    rsi = np.zeros_like(prices)
    rsi[:n] = 1. - 1./(1. + rs)

    for i in range(n, len(prices)):
        delta = deltas[i - 1]  
        
        if delta > 0:
            upval = delta
            downval = 0.
        else:
            upval = 0.
            downval = -delta

        up = (up*(n - 1) + upval)/n
        down = (down*(n - 1) + downval)/n

        rs = up/down
        rsi[i] = 1. - 1./(1. + rs)

        data_merged['RSI'] = rsi

relative_strength(data_merged['SPY_prior'], n=14)

#MACD indicator
def moving_average(x, n=9):
    sma = pd.rolling_mean(x, n)
    return sma

def moving_average_convergence(x, nslow=26, nfast=12):
    emaslow = pd.ewma(x, span=nslow, min_periods=1)
    emafast = pd.ewma(x, span=nfast, min_periods=1)
    result = pd.DataFrame({'MACD': emafast-emaslow, 'emaSlw': emaslow, 'emaFst': emafast})
    return result

data_merged[['MACD', 'emaSlw', 'emaFst']] = moving_average_convergence(data_merged['SPY_prior'], 26, 12)
data_merged['ema_9_MACD'] = pd.ewma(data_merged['MACD'], span=9, min_periods=1)
#when the MACD crosses above the 9 ema of the MACD, it is considered bullish and the opposite is bearish. let's take the difference between MACD and ema_9_MACD to show this
data_merged['MACD_signal'] = data_merged['MACD'] - data_merged['ema_9_MACD']

#ROC indicator
def ROC(data_merged, n):  
    M = data_merged['SPY_prior'] - data_merged['SPY_prior'].shift(n)
    N = data_merged['SPY_prior'].shift(n)  
    ROC = pd.Series(M / N) 
    data_merged["ROC"] = ROC 
    
#Calculating ROC    
ROC(data_merged, 14)

#Stochastic oscillator.
#Stochastic Oscillator is a momentum indicator that shows the location of the close(price at the end of the day) relative to the high-low range over a set number of periods.
 
#creating the function.
def STO(data_merged, n):  
    SOk = (data_merged['SPY_prior'] - data_merged['Low_prior'].rolling(window=n).min()) / (data_merged[
                'High_prior'].rolling(window=n).max() - data_merged['Low_prior'].rolling(window=n).min()) 
    
    SOd = pd.ewma(SOk, span = 3) 
    data_merged['S0d'] = SOd  
    
#Calculating STO
STO(data_merged, 14)


#Internal Bar Strength Indicator. 
#It takes a value of 0 if the closing (price at the end of the day) price is the lowest price of the day, and 1 if the closing price is the highest price of the day.
#I found this blog and it does a pretty good job in predicting the market - http://jonathankinlay.com/2016/06/the-internal-bar-strength-indicator/
#IBS  =  (Close – Low) / (High – Low)
#there is actually a short paper on this - http://qusma.com/wp-content/uploads/2013/09/The-IBS-Effect-Mean-Reversion-in-Equity-ETFs.pdf

data_merged['IBS'] = (data_merged['SPY_prior'] - 
                      data_merged['Low_prior']) / (data_merged['High_prior'] - data_merged['Low_prior'])

#Let's work with the volume data
#Excess volume is the difference between the daily volume and the 30-day moving average volume expressed as a percentage of the 30-day moving average. 
#If this number is positive, the daily volume for the given day is higher than the 30-day moving average. Otherwise, it is below. 
data_merged['excess_volume_prior'] = (data_merged['Volume_prior'] - pd.rolling_mean(data_merged['Volume_prior'], window = 30))/ (pd.rolling_mean(data_merged[
            'Volume_prior'], window = 30))

#We have one more price based indicator to calculate. That is, Z-score. 
#In order to prevent data snooping, we need to make sure we do expanding z-score. 
#This means that we use the data we have at a given point in time to calculate the mean and the standard deviation.

data_merged['SPY_Z_score'] = (data_merged['SPY_prior_ret'] - data_merged[
        'SPY_prior_ret'].expanding(min_periods=1).mean())/data_merged['SPY_prior_ret'].expanding(min_periods=1).std()

#let's make another column that shows the direction of the daily market movement. 
#If the percentange change is negative (this means that the price is lower), we will have -1, if it is positive, we will have 1. 0 return will be given -1. 
#This will help in building classification models. 

#define target direction
data_merged['SPY_direction'] = np.where(data_merged['SPY_ret'] > 0, 1, -1)
data_merged['SPY_direction'] = data_merged['SPY_direction'].astype('category')

#Let's get rid of rows that have NAs in it. 
data_merged.dropna(inplace=True)

#let's apply a non-visual way of checking for NAs
print(np.any(np.isnan(data_merged.ix[:, 2:32])))
#I get FALSE for this code and it shows that we don’t have any NA value in any of the rows

#let's check if we have any infinite values. 
print(np.all(np.isfinite(data_merged.ix[:,2:32])))
#I get TRUE for this code indicating that all the values are finite

#Along the way, we have dropped some rows. So, our dataframe indexing isn’t starting from 0. Let’s fix this. This will become handy when we are doing time series split.
data_merged.drop('index', axis = 1, inplace = True)
data_merged.reset_index(inplace=True)
data_merged.drop('index', axis=1,inplace=True)

#let's see how many observations we have. 
len(data_merged)

#Let's create training and testing subsets
cutoff = 4000

train_set = pd.DataFrame()
test_set = pd.DataFrame()
train_set = data_merged[data_merged.index < cutoff].copy()
test_set = data_merged[data_merged.index >= cutoff].copy()

print('Observations: %d' % (len(data_merged)))
print('Training Observations: %d' % (len(train_set)))
print('Testing Observations: %d' % (len(test_set)))
plt.plot(train_set['SPY_ret'])
plt.plot([None for i in train_set['SPY_ret']] + [x for x in test_set['SPY_ret']])
plt.show()
#the blue ones are in the training set and the green ones in the testing set. 

#Let's seperate the features and the target. 
# training features
train_features = train_set[['SPY_prior_ret','VIX_prior_ret','SPY_trend','RSI','MACD_signal', 'ROC','S0d',
                            'SPY_Z_score','excess_volume_prior','IBS','MSCI_EM_prior_ret','DP_prior', 'PE_prior', 'BM_prior',
                            'CAPE_prior','BY_prior', 'DEF_prior','TERM_prior', 'CAY_prior', 'SIM_prior', 
                            'VRP_prior', 'IC_prior', 'BDI_prior','NOS_prior', 'CPI_prior', 'PCR_prior', 'MA_prior',
                            'OIL_prior', 'SI_prior']]

# training target 
train_target = train_set['SPY_ret']

#testing features
test_features = test_set[['SPY_prior_ret','VIX_prior_ret','SPY_trend','RSI','MACD_signal', 'ROC','S0d','SPY_Z_score',
                          'excess_volume_prior','IBS','MSCI_EM_prior_ret','DP_prior', 'PE_prior', 
                          'BM_prior', 'CAPE_prior','BY_prior', 'DEF_prior','TERM_prior', 'CAY_prior', 'SIM_prior',
                          'VRP_prior', 'IC_prior', 'BDI_prior','NOS_prior', 'CPI_prior', 'PCR_prior', 'MA_prior', 
                          'OIL_prior', 'SI_prior']]
#'PCA.tech_prior'

#testing target 
test_target = test_set['SPY_ret']

#pick out the dates for the test set
dates = test_set['Period']

def fitAndPlot(train_features, test_features, model, **extraArgs):

    # Create linear regression object
    SPY_ret_model = model(**extraArgs)

    # Fit the data
    SPY_ret_model.fit(train_features,train_target)

    # Get the predictions
    predictions = SPY_ret_model.predict(test_features)
    
    #let's design the strategy
    directions = pd.DataFrame()
    directions = np.sign(predictions)
    strategy = directions * test_target
    
    strategy_equity_curve = strategy.cumsum().apply(np.exp)
    buy_and_hold_equity_curve = test_target.cumsum().apply(np.exp)

    # Plot the actuals and the predictions
    f, (ax1) = plt.subplots(1, 1, figsize=(7,5))

    ax1.plot(dates,strategy_equity_curve,linestyle='--', label = 'Strategy')
    ax1.plot(dates,buy_and_hold_equity_curve, label = 'SPY_buy&hold')
    ax1.legend(loc=0)
    ax1.set_ylabel('Performance')
    ax1.set_xlabel('Dates')
    f.autofmt_xdate()
    #plt.tight_layout()    
 

    print("Model RMS Error: {0:.10f} for {1}".format(np.sqrt(np.mean((predictions - test_target) ** 2)),model.__name__))
    print('Annualized return') 
    print((strategy.mean()*252)*100)
    print('Volatility') 
    print((strategy.std()*252**0.5)*100)
    print('Sharpe Ratio')
    print((strategy.mean()*252) / (strategy.std()*252**0.5))
    
    return (SPY_ret_model)

fitAndPlot(train_features, test_features, LinearRegression)

# Linear Regression do well when the features are uncorrelated or have low correlation. So, let's try to filter the features by 
# looking at the correlation. 

correl = data_merged[['SPY_prior_ret','VIX_prior_ret','SPY_trend','RSI','MACD_signal', 'ROC','S0d','SPY_Z_score','excess_volume_prior','IBS',
                           'MSCI_EM_prior_ret','DP_prior', 'PE_prior', 'BM_prior', 'CAPE_prior', 'PCA.price_prior', 'BY_prior', 'DEF_prior',
                            'TERM_prior', 'CAY_prior', 'SIM_prior', 'VRP_prior', 'IC_prior', 'BDI_prior','NOS_prior', 
                            'CPI_prior', 'PCR_prior', 'MA_prior','OIL_prior', 'SI_prior']].corr()

# Generate a mask for the upper triangle
mask = np.zeros_like(correl, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True

# Set up the matplotlib figure
f, ax = plt.subplots(figsize=(20, 20))

# Generate a custom diverging colormap
cmap = sns.diverging_palette(220, 10, as_cmap=True)

# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(correl, annot=True, mask=mask, cmap=cmap, vmax=.3,
            square=True,linewidths=.5, cbar_kws={"shrink": .5}, ax=ax)

#SPY_prior_ret has higher correlation with the technical indicators than with the fundamental indicators. 
#Let's see how the linear model does when we use only the prior return and the fundamental data. 
#training features
new_train_features = train_set[['SPY_prior_ret','DP_prior', 'PE_prior', 'BM_prior', 'CAPE_prior', 
                                'BY_prior', 'DEF_prior','TERM_prior', 'CAY_prior', 
                                'BDI_prior','NOS_prior','CPI_prior','SI_prior']]

#testing features
new_test_features = test_set[['SPY_prior_ret','DP_prior', 'PE_prior', 'BM_prior', 'CAPE_prior',
                              'BY_prior', 'DEF_prior','TERM_prior', 'CAY_prior',
                              'BDI_prior','NOS_prior','CPI_prior','SI_prior']]

 fitAndPlot(new_train_features,new_test_features,LinearRegression)

from sklearn import feature_selection
from sklearn.feature_selection import SelectKBest
from sklearn.feature_selection import chi2

select_feature = feature_selection.SelectKBest(score_func=feature_selection.f_regression,k=10)
results = select_feature.fit(train_features, train_target)
print(results.scores_)
print(results.pvalues_)

#training features
newer_train_features = train_set[['SPY_prior_ret', 'SPY_Z_score', 'MSCI_EM_prior_ret', 'IBS', 'MACD_signal']]

#testing features
newer_test_features = test_set[['SPY_prior_ret', 'SPY_Z_score', 'MSCI_EM_prior_ret', 'IBS', 'MACD_signal']]

fitAndPlot(newer_train_features, newer_test_features, LinearRegression)

all_train_features = train_set[['SPY_prior_ret','VIX_prior_ret','SPY_trend','RSI','MACD_signal', 'ROC','S0d',
                            'SPY_Z_score','excess_volume_prior','IBS','MSCI_EM_prior_ret','DP_prior', 'PE_prior', 'BM_prior',
                            'CAPE_prior','BY_prior', 'DEF_prior','TERM_prior', 'CAY_prior', 'SIM_prior', 
                            'VRP_prior', 'IC_prior', 'BDI_prior','NOS_prior', 'CPI_prior', 'PCR_prior', 'MA_prior',
                            'OIL_prior', 'SI_prior','PCA.tech_prior','PCA.price_prior']]

#testing features
all_test_features = test_set[['SPY_prior_ret','VIX_prior_ret','SPY_trend','RSI','MACD_signal', 'ROC','S0d','SPY_Z_score',
                          'excess_volume_prior','IBS','MSCI_EM_prior_ret','DP_prior', 'PE_prior', 
                          'BM_prior', 'CAPE_prior','BY_prior', 'DEF_prior','TERM_prior', 'CAY_prior', 'SIM_prior',
                          'VRP_prior', 'IC_prior', 'BDI_prior','NOS_prior', 'CPI_prior', 'PCR_prior', 'MA_prior', 
                          'OIL_prior', 'SI_prior','PCA.tech_prior','PCA.price_prior']]

#After running Random Forest Regression, AdaBoost Regression, and Gradient Boosting Regression models, the only model that beats the Benchmark Linear Model is the Random Forest model. 
a = np.array(['Random Forest Regression','AdaBoost Regression', 'Gradient Boosting'], dtype=np.str_)
b = ['17.52%','0.39%','5.75%']
c = ['15.14%','15.18%','15.17%']
d = ['1.16','0.025','0.38']
t = Table([a, b, c,d], names=('Model type ', 'Annualized Return', 'Annualized Volatility ', 'Sharpe Ratio '))
t

#Let's first try to optimize random forest regression. 
min_leaf_grid = range(1,80,4)
rms_scores = []
for  min_leaf in min_leaf_grid:
    
    print("Working min_leaf {}".format(min_leaf))
    SPY_ret_model = RandomForestRegressor(min_samples_leaf=min_leaf, n_estimators=300,random_state=32)

    # Fit the data
    SPY_ret_model.fit(all_train_features,train_target)

    # Get the predictions
    predictions = SPY_ret_model.predict(all_test_features)
    
    rms_score = np.sqrt(np.mean((predictions - test_target) ** 2))
    rms_scores.append(rms_score)

plt.plot(min_leaf_grid, rms_scores,marker='o')
plt.xlabel("Minimum Number of Leaves")
plt.ylabel("RMS Error")

fitAndPlot(all_train_features,all_test_features, RandomForestRegressor, n_estimators=1000, min_samples_leaf=30,random_state=32)

a = np.array(['Optimized Random Forest Regression'], dtype=np.str_)
b = ['18.73%']
c = ['15.14%']
d = ['1.24']
t = Table([a, b, c,d], names=('Model type ', 'Annualized Return', 'Annualized Volatility ', 'Sharpe Ratio '))
t