<jupyter_start><jupyter_text>Loading up required libraries and configurations<jupyter_code>import quandl
import pandas_datareader.data as web
import datetime
import pandas as pd
import sklearn
import numpy as np
from collections import defaultdict
from IPython.display import display
import scipy as sp
from operator import methodcaller

# evaluate usage
pd.options.mode.chained_assignment = None  # default='warn'
""""
usage of this API key is monitored
please don't use this key for any other work, neither make it available on the web by any means
if you would like to access the same API for a different project,
please create an account in quandl.com (it is free) and generate your own API key
""" 
quandl.ApiConfig.api_key = "<KEY>"<jupyter_output><empty_output><jupyter_text>Getting the data  Stock data For the stock market general data, we will use Yahoo API, which contains reliable free access data for a number of stock markets, including Bovespa.<jupyter_code>start_date = '1998-1-1'
end_date = '2014-12-31'
symbol = "PETR4.SA"
df = web.DataReader(symbol, 'yahoo', start_date, end_date)<jupyter_output><empty_output><jupyter_text>General Market Data For the general market data, there is excelent data made available by Banco Central do Brasil (brazilian central bank). The data is available through a webservice, but there is a neat API quandl which makes it easier to access this same data. With a free profile we have limited access, but it is enough for the following tests.There are over 1,300 different indicators available, from different time periods including daily, weekly, monthly, yearly, etc. At the moment we will stick with 10 relevant indicators which are available daily, and then move on to add more data as we see fit.<jupyter_code>daily = {
    'Selic': 'BCB/432',
    'Exchange Rate USD Sell': 'BCB/1',
    'Exchange Rate USD Buy': 'BCB/10813',
    'BM&F Gold gramme': 'BCB/4',
    'Bovespa total volume': 'BCB/8',
    'International Reserves': 'BCB/13621',
    'Bovespa index': 'BCB/7',
    'Foreign exchange operations balance': 'BCB/13961',
    'Nasdaq index': 'BCB/7810',
    'Dow Jones index': 'BCB/7809'
}

## removed montly indicators for now - to be added later
# monthly = {
#     'IPCA-E': 'BCB/10764',
#     'IGP-M': 'BCB/7451',
#     'IPCA-15': 'BCB/74/78',
#     'Net public debt': 'BCB/4473'
# }
for var, code in daily.items():
    df = pd.concat([df, quandl.get(code, start_date=start_date , end_date=end_date)], join='inner', axis=1)
    df = df.rename(columns={'Value': var})<jupyter_output><empty_output><jupyter_text>Trend indicators The trend indicators are borrowed from the field known as technical analysis, or graphism, that aims to find patterns by analyzing the trend of price and volume.We will start with the most known: Moving Averages (indicator: MACD, moving average convergence divergence), Momentum, Daily Returns,  Already included in the dataset* Momentum: momentum is nothing but the current price, divided by the price X days earlier. The momentum is already included the dataset when we analyse the trend for Adj Close and all other variables* Daily Return: it is the same as the momentum, but for one day before. To include:* Moving Average: moving average for X days. important to understand longer term trends* Bollinger Bands: 95% confidence interval for the moving averages.* CandleStick* Elliot Waves* Volume/Price<jupyter_code># moving average and bollinger bands
for n in range(10,61,10):
    df['sma'+str(n)] = df['Adj Close'].rolling(window=n, center=False).mean()
    std =df['Adj Close'].rolling(window=n, center=False).std()
    df['bb_lower'+str(n)] = df['sma'+str(n)] - (2 * std)
    df['bb_upper'+str(n)] = df['sma'+str(n)] + (2 * std)<jupyter_output><empty_output><jupyter_text>Creating the labels The general approach to the stock market problem is to use non-linear regressors to predict future prices. Although it is easy to predict the price for the day ahead, as you move days further the r2_score sinks and the prediction becomes useless. The inovative approach we will follow here is meant to treat this problem as a classification problem. In order to treat this problem as a classifier, we will pre-define a trading strategy, backtest it to the full dataset, and define the label to be whether this trading strategy results in a successful trade or not.  Swing Trade The first and most simple analysis is a swing trade strategy. We will buy the stock, and hold it for n days. If it reaches an upper boundary (+x%), we sell. If it reaches a lower boundary (-y%), we will also short our position.So the challenge is within the next n days the stock needs to reach the upper boundary, before it reaches the lower boundary. That will be counted as a successful trade. If it reaches the lower boundary before, or n days passes without reaching the upper boundary, the trading is considered a failure.The name swing trade means we are speculating on the week/bi-week pattern, the swing, not necessarily on the longer term trend.The parameters n, x, and y, will be optimized through a genetic algorithm to achieve the optimal trading strategy for this setup.<jupyter_code>def create_labels(input_df, forward=10, profit_margin=.05, stop_loss=.05):
    
    df = input_df.copy()
    for row in range(df.shape[0]-forward):

        # initialize max and min ticks
        max_uptick = 0
        min_downtick = 0 

        # move all days forward
        for i in range(1,forward+1):
            delta = (df.ix[row+i, 'Adj Close'] / df.ix[row, 'Adj Close'])-1
            if delta > max_uptick:
                max_uptick = delta
            if delta < min_downtick:
                min_downtick = delta

        # evaluate ticks against predefined strategy parameters
        if max_uptick >= profit_margin and min_downtick <= -stop_loss:
            df.ix[row,'Label'] = 1
        else:
            df.ix[row,'Label'] = 0        

    return df.dropna()

df = create_labels(df)<jupyter_output><empty_output><jupyter_text>Rounding up the features Now we've got the main data, and the label, let's start creating our features. There is also some innovation going on here, we are running from the traditional approach of stock price prediction. The first and foremost difference is that instead of analyzing the raw data, I want to analyze the trend for each variable. So in a single day I will look at how the variable increased or decreased in the last N days. To center at 0, for each variable I will calculate (value / value in n-1). The number of days I will look back may also vary, and it will also depend on the trading strategy to follow. But for simplicity we will start with a round number such as 60 days for the swing trading strategy of around 10 days.<jupyter_code>def create_features(input_df, base = 60):
    """ Receives a dataframe as input 
        Returns a new dataframe with ratios calculated
    """
    # avoid modifying in place
    df = input_df.copy() 
    # get all columns ubt the label
    cols = list(df.columns)
    cols.remove('Label')
    # create the new columns for the number of days
    for n in range(1,base+1):
        new_cols = list(map(lambda x: "{}-{}".format(x,n), cols))
        df[new_cols] = (df.loc[:, cols] / df.shift(n).loc[:, cols]) - 1
    # remove inf and drop na
    df.replace([np.inf, -np.inf], np.nan, inplace=True) 
    df.dropna(axis=0, inplace=True)
   
    # leave or remove original columns? for now I will leave them
    #return df.drop(cols, axis=1)
    return df
df = create_features(df)
df.shape<jupyter_output><empty_output><jupyter_text>Understanding the label and features Let's start by analizying the label distribution. This will tell us a lot about the dataset, which optimal classifier to choose, and whether we will need to use a stratified approach when splitting the dataset for testing or cross-validation.<jupyter_code># break up X and y arrays, convert to numpy array
features = [x for x in df.columns if x != "Label"]
X = df[features].values
y = df['Label'].values
X.shape, y.shape
np.bincount(y.astype(int)) / len(y)<jupyter_output><empty_output><jupyter_text>As expected, there will be only a few occurrences where such a trading strategy results in success. Only 6,25% of the observations have label 1(success), while 93,75% have label 0 (failure). A stratified approach will be required when splitting the dataset later. Next step is to take a look at the features we have. I will start standardizing to z-score, then checking and removing outliers, and finally analyse which features are most relevant and attempt to extract principal components. For now, I will keep working with the full data, since I'm looking for understanding data, I'm doing predictions at this point. I will come back later to this point to divide the dataset<jupyter_code># scaling
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X = scaler.fit_transform(X)
# feature selection
from sklearn.feature_selection import SelectKBest
f_selector = SelectKBest()
f_selector.fit(X, y)
sorted(zip(features, f_selector.scores_, f_selector.pvalues_), key=lambda x: -x[1])<jupyter_output><empty_output><jupyter_text>Strangely enough, it shows an incredible high importance for International Reserves features, as well for other features related to the outside market, such as Nasdaq, and Exchange Rate.Petrobras, the company being analyzed, is very influenced by the price and availability of USD in brazilian market. It is debts are most in dollar and the price of its main product, oil, is set globally. So these findings that my be accurate, just random noise, or due to some mistake in the data preprocessingWill come back later to this point if this correlation is not confirmed with other analysis.<jupyter_code># pca conversion
# okay, I've got all these data points
from sklearn.decomposition import PCA
pca = PCA(n_components = 32)
pca.fit(X)
sum(pca.explained_variance_ratio_)

# why 32? I'm aiming at 80% explained variance, with the minimum number of components
# I just tried a few and got to the required number of components
# this could also be optimized later, aimed at performance
X_transformed = pca.transform(X)<jupyter_output><empty_output><jupyter_text>Okay, now we've transformed into fewer principal components, thus helping to avoid the curse of dimensionality, let's take a step back and see which are the main variables impacting each of the first components. Are they the same variables which stand out in the Anova F-value tests conducted before?<jupyter_code># checking the variables with most impact on the first component
i = np.identity(len(features)) # identity matrix
coef = pca.transform(i)
sorted(zip(features, coef[:, 1]), key=lambda x:-x[1])<jupyter_output><empty_output><jupyter_text>The results hardly show something relevant. The coefficientes between the features are very similar for the first principal component<jupyter_code># I will try feature selection again, with the components
from sklearn.feature_selection import f_classif
f_selector = SelectKBest(f_classif)
f_selector.fit(X_transformed,y)
sorted(zip(f_selector.scores_, f_selector.pvalues_), key=lambda x: -x[0])<jupyter_output><empty_output><jupyter_text>The p values for the Anova F-value test, all above 0.5, tell a worrisome story. None has a verifiable correlation with the label. Nevertheless, let's push forward and see how good a prediction we can make with this data  Predicting Let's move on to prediction and see which results we can get with the work that is already done.A previous work was done using several classifiers, and Gradient Boosting was the most promising. It is a model that can fit well complex models, and is one of the top perfomer algorithms used in the wide. Let's start with it and check the results<jupyter_code>from sklearn.model_selection import StratifiedShuffleSplit
from sklearn.ensemble import GradientBoostingClassifier as GBC
from sklearn.neighbors import KNeighborsClassifier as kNN
from sklearn.model_selection import cross_val_score
import warnings; warnings.filterwarnings("ignore")
# define cross validation strategy
cv = StratifiedShuffleSplit(n_splits=10, test_size=.1, random_state=42)
# initiate, train and evaluate classifier
clf = GBC(random_state=42)
scores = cross_val_score(clf, X, y, cv=cv, scoring='precision')
print("Precision: %0.2f (+/- %0.2f)" % (scores.mean(), scores.std() * 2))
# same, with kNN
clf = kNN()
scores = cross_val_score(clf, X, y, cv=cv, scoring='precision')
print("Precision: %0.2f (+/- %0.2f)" % (scores.mean(), scores.std() * 2))<jupyter_output>Precision: 0.17 (+/- 0.63)<jupyter_text>40% is a reasonable precision, considering it is a highly skewed dataset, with only 6,25% of the observations with label 'Succeess'. But still a lot of variation - with a 95% confidence interval, the precision can be as low as .17. Let's check with the principal components, to see if we can get better results.I will try with different number of components<jupyter_code># components = [20,50,100,150,200,300,400,len(features)]
components = [10,20,50,100, 200, 300]
for n in components:
    # pca
    pca = PCA(n_components = n, random_state=42)
    X_transformed = pca.fit_transform(X)
    
    # predict
    clf = GBC(random_state=42)
    scores = cross_val_score(clf, X_transformed, y, cv=cv, scoring='precision')
    print("Precision {} components: {:0.2f} (+/- {:0.2f})".format(
            n, scores.mean(), scores.std() * 2))
# same, with knn
components = [10, 20, 50, 100, 200, 300]
for n in components:
    # pca
    pca = PCA(n_components = n, random_state=42)
    X_transformed = pca.fit_transform(X)
    
    # predict
    clf = kNN()
    scores = cross_val_score(clf, X_transformed, y, cv=cv, scoring='precision')
    print("Precision {} components: {:0.2f} (+/- {:0.2f})".format(
            n, scores.mean(), scores.std() * 2))<jupyter_output>Precision 10 components: 0.05 (+/- 0.30)
Precision 20 components: 0.07 (+/- 0.32)
Precision 50 components: 0.07 (+/- 0.32)
Precision 100 components: 0.07 (+/- 0.32)
Precision 200 components: 0.12 (+/- 0.60)
Precision 300 components: 0.23 (+/- 0.79)<jupyter_text>The results with the principal components doesn't seem any better. I've also tried with 20 folds in the cross validation to reduce the margin of error, and between 100 to 200 components I get a precision of around .38 and a margin of error of around .5, which too high.  Optimizing parameters with genetic algorithm Well, next step is to change my strategy. In the swing trade strategy, if I set number of days 10, profit margin(x) to 5% and stop loss(y) to 5%, I have to be right at least 51% of the times to actually win some money. That is not the case, as my current precision lower boundary is 17%.But I might have a better chance identifying a different variation of this strategy. So as discussed before, let's try to optimize these 3 parameters: n, x and y. We will use a genetic algorithm which will use precision score from the cross validation function as its own scoring function, to determine which variation will perpetuate.  In order to that, we will organize the code to create the dataset into classes and modules, along with the Genetic Algoritm, in the attached files stock.py and strategies.py, and import them to the notebook.<jupyter_code>import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning)

class Backtester():

    def __init__(self, df):
        self.df = df

    def create_labels(self, forward=10, profit_margin=.05, stop_loss=.05):

        for row in range(self.df.shape[0]-forward):

            # initialize max and min ticks
            max_uptick = 0
            min_downtick = 0 

            # move all days forward
            for i in range(1,forward+1):
                delta = (self.df.ix[row+i, 'Adj Close'] / self.df.ix[row, 'Adj Close'])-1
                if delta > max_uptick:
                    max_uptick = delta
                if delta < min_downtick:
                    min_downtick = delta

            # evaluate ticks against predefined strategy parameters
            if max_uptick >= profit_margin and min_downtick <= -stop_loss:
                self.df.ix[row,'Label'] = 1
            else:
                self.df.ix[row,'Label'] = 0        

        self.df.dropna(inplace=True)

    def prep_data(self):
        features = [x for x in self.df.columns if x != "Label"]
        X = self.df[features].values
        y = self.df['Label'].values
        return X,y

    def score(self, X, y):
        # apply PCA
        pca = PCA(n_components = 10, random_state=42)
        X_transformed = pca.fit_transform(X)
        
        #predict
        clf = GBC(random_state=42)
        cv = StratifiedShuffleSplit(n_splits=10, test_size=.1, random_state=42)
        scores = cross_val_score(clf, X_transformed, y, cv=cv, scoring='precision')

        # return score
        return (scores.mean())

    def evaluate(self, forward, profit_margin, stop_loss):
        self.create_labels(forward=forward, profit_margin=profit_margin, stop_loss=stop_loss)		
        score = self.score(*self.prep_data())
        print("span: {}, profit_margin: {:.3f}, stop_loss: {:.3f} --  score: {:.3f}".format(
            forward, profit_margin, stop_loss, score))
        return score
class Strategy():

    def __init__(self, span = 7, profit_margin = .06, stop_loss = .04):
        self.span = span
        self.profit_margin = profit_margin
        self.stop_loss = stop_loss
        self.mutations = [
            self.increase_span,
            self.decrease_span,
            self.increase_stop_loss,
            self.decrease_stop_loss,
            self.increase_profit_margin,
            self.decrease_profit_margin
        ]

    def mutate(self):
        np.random.choice(self.mutations)()

    def inform_params(self):
        return self.span, self.profit_margin, self.stop_loss

    def report(self):
        print("span: {}, profit_margin: {:.3f}, stop_loss: {:.3f}".format(
            self.span, self.profit_margin, self.stop_loss))

    # add a random component to mutation
    # allow 'wild' mutations

    def increase_span(self):
        self.span += 2

    def decrease_span(self):
        self.span -= 2

    def increase_profit_margin(self):
        self.profit_margin = self.profit_margin * 1.3

    def decrease_profit_margin(self):
        self.profit_margin = self.profit_margin * .7

    def increase_stop_loss(self):
        self.stop_loss = self.stop_loss * 1.3

    def decrease_stop_loss(self):
        self.stop_loss = self.stop_loss * .7


class GA():

    def __init__(self, df):
        """ Seed 2 initial strategies and an instance of backtester """

        self.backtester = Backtester(df.copy())

        self.strategies = pd.DataFrame(np.zeros((2,2)), columns = ['strat', 'score'])
        self.strategies['strat'] = self.strategies['strat'].apply(lambda x:Strategy())
        self.strategies['score'] = self.strategies['strat'].apply(self.evaluate)


    def fit(self, cycles):
        """ Run evolution for n cycles """
        i = 0
        while i < cycles:
            self.reproduce()
            # self.select()

            i += 1

    def best_strategy(self):
        """ Sort and return top perform in available strategies """

        self.strategies =  self.strategies.sort_values(by='score', ascending=False)
        self.strategies.iloc[0, 0].report()
        print("score: {:.4f}".format(self.strategies.iloc[0, 1]))

    def evaluate(self, strat):
        """ To implement: 
            Should evaluate only for those which value is zero, to avoid the cost of re-evaluating 
        """
        return self.backtester.evaluate(*strat.inform_params())

    def reproduce(self):
        """ Create new strategy based on its parents. """

        # sort and take top two performers in the list
        parents = self.strategies.sort_values(by='score', ascending=False).iloc[:2, 0]

        # create six offsprings
        for _ in range(6):
            stratN = self.crossover(*parents)
            stratN.mutate()
            # setting with enlargement using index based selection (not available for position based)
            self.strategies.ix[self.strategies.shape[0]] = (stratN, self.evaluate(stratN))

            # remove identical offspring, there is no use

    def crossover(self, stratA, stratB):
        """ Choose between parents attributes randomly. Return new strategy """

        span = np.random.choice([stratA.span, stratB.span])
        stop_loss = np.random.choice([stratA.stop_loss, stratB.stop_loss])
        profit_margin = np.random.choice([stratA.profit_margin, stratB.profit_margin])
        return Strategy(span=span, stop_loss=stop_loss, profit_margin=profit_margin)

    def select(self):
        """ Remove strategies which are bad performers 
            Define bad as 50 percent worst than best """

        # define cut off score as 50% of the max score
        cut_off = self.strategies['score'].max() * .75
        # remove strategies with scores below the cut-off
        self.strategies = self.strategies[self.strategies['score'] >= cut_off]


# ga = GA(df)
# ga.fit(20)
# ga.best_strategy()<jupyter_output>span: 7, profit_margin: 0.060, stop_loss: 0.040 --  score: 0.000
span: 7, profit_margin: 0.060, stop_loss: 0.040 --  score: 0.000
span: 7, profit_margin: 0.060, stop_loss: 0.052 --  score: 0.000
span: 9, profit_margin: 0.060, stop_loss: 0.040 --  score: 0.100
span: 7, profit_margin: 0.060, stop_loss: 0.028 --  score: 0.000
span: 5, profit_margin: 0.060, stop_loss: 0.040 --  score: 0.000
span: 7, profit_margin: 0.060, stop_loss: 0.052 --  score: 0.000
span: 9, profit_margin: 0.060, stop_loss: 0.040 --  score: 0.100
span: 11, profit_margin: 0.060, stop_loss: 0.040 --  score: 0.000
span: 11, profit_margin: 0.060, stop_loss: 0.040 --  score: 0.000
span: 9, profit_margin: 0.060, stop_loss: 0.052 --  score: 0.000
span: 9, profit_margin: 0.042, stop_loss: 0.040 --  score: 0.153
span: 11, profit_margin: 0.060, stop_loss: 0.040 --  score: 0.000
span: 9, profit_margin: 0.060, stop_loss: 0.052 --  score: 0.000
span: 9, profit_margin: 0.029, stop_loss: 0.040 --  score: 0.025
span: 9, profit_margin[...]<jupyter_text>GA is not helping too much. I will leave this as an option to optimize in the end, but for now, I will go back to the basics.I went in the PCA line. Let's try feature selecting and see what results we can get  Feature Selection How many features is the optimal number?Which classifier should I use?Which parameter can I get?<jupyter_code># results with 10 folds (stratified), precision<jupyter_output><empty_output><jupyter_text>Removed code<jupyter_code># this is just a test.
# it shows the proportions have been kept in the training label
# gotta find another explanation for the errors show in the PCA
cv.get_n_splits(X_transformed, y)
for train_index, test_index in cv.split(X_transformed, y):
    y_train, y_test = y[train_index], y[test_index]
    print(np.bincount(y_train.astype(int)))
    print(np.bincount(y_test.astype(int)))
# outlier detection
# vectorized
# calculate step = (q3-q1*1.5)

df.quantile(.5)
(df.quantile(.75) - df.quantile(.25))*5

# don't remove outliers for now. come back to this point later
# okay, I've got a all this data - is a lot - I need to drop the original columns, they are irrelevant. 
# drop or scale. I will drop for now
for a,b in new_row.items():
    print(a)
#testing
#inspection_cols = ['Open', 'High', 'Low']
#cols = list(inspection_cols) + \
#        list(map(lambda x: "{}-{}".format(x,1), inspection_cols)) + \
#        list(map(lambda x: "{}-{}".format(x,2), inspection_cols))
# will run forever. this will not work
def create_features(df, base = 7):
    """ - Get a basis, X days before, and calculate all the other variables on that basis
        - Lookup_idx calculation: will look to 7 days before, including this one.
            7 - (7-1) = 7 - 6 = 1
            7 - (7-2) = 7 - 5 = 2
            ...
            7 - (7-7) = 7 -0 = 7  
            
        - Should not return n/as
    """
    for idx, row in enumerate(df.iterrows()):
        if idx>=base:
            base_row = df.iloc[idx-base, :]
            for n in range(1,base+1):
                lookup_idx = idx-(base-n)
                new_row = df.iloc[lookup_idx, :] / base_row
                for col, value in new_row.items():
                    df.ix[idx, col+"-"+str(n)] = value
    return df

# check and remove outliers
# use Turkey's method - 1,5 * IQR

outliers_count = defaultdict(int) 
for feature in df.columns:    
    Q1 = np.percentile(df[feature], 25)
    Q3 = np.percentile(df[feature], 75)
    step = (Q3-Q1)*1.5
    
    # Display the outliers
    print("Data points considered outliers for the feature '{}':".format(feature))
    outliers_for_feature = df[~((df[feature] >= Q1 - step) & (df[feature] <= Q3 + step))]
    for idx in  outliers_for_feature.index.values:
        outliers_count[idx] += 1
    
    (outliers_for_feature.shape)
    
# print the total list of outliers
print("Total number of data points identified as outliers: {}".format(len(outliers_count)))

outliers_count
# OPTIONAL: Select the indices for data points you wish to remove
#outliers =  list(map(lambda (x,y):x, filter(lambda (x,y):y>1, outliers_count.items())))
#print("Total number of data points identified as outliers in more than one data entry: {}".format(len(outliers)))

# Remove the outliers, if any were specified
# good_data = log_data.drop(log_data.index[outliers]).reset_index(drop = True)
#print(outliers)<jupyter_output>Data points considered outliers for the feature 'Open-1':
{'fraction_non_alphanumeric': 0.09795281731331644, 'fraction_numerical': 0.02398128290115032, 'mean_word_length': 3.548776161759489, 'pattern_counts': {'":': 0, '<': 68, '<?xml version=': 0, '>': 70, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '200427', 'license': 'MIT', 'n_tokens_mistral': 8247, 'n_tokens_neox': 7521, 'n_words': 2909, 'path': 'ENIAC-experiments/v-previous/quandl-v2.ipynb', 'repo_name': 'lucasosouza/oraclia-research', 'size': '25660'}