%config InlineBackend.figure_format = 'retina'

import numpy as np
import pandas as pd
from scipy.special import logit
from datetime import datetime
from dateutil.relativedelta import relativedelta
from tech import technical
import warnings
warnings.filterwarnings('ignore')

import matplotlib.pyplot as plt
import seaborn as sns

from xgboost import XGBClassifier
from sklearn.decomposition import PCA, KernelPCA
from sklearn.cross_validation import KFold, cross_val_score
from sklearn.metrics import make_scorer, accuracy_score, precision_score, recall_score
from sklearn.grid_search import GridSearchCV
from sklearn.feature_selection import VarianceThreshold, RFE, SelectKBest, chi2
from sklearn.preprocessing import MinMaxScaler
from sklearn.pipeline import Pipeline, FeatureUnion
from sklearn.linear_model import LogisticRegression
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.neighbors import KNeighborsClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.svm import SVC
from sklearn.ensemble import BaggingClassifier, ExtraTreesClassifier, GradientBoostingClassifier, VotingClassifier, RandomForestClassifier, AdaBoostClassifier

sns.set_style('whitegrid')
pd.set_option('display.max_columns', None) # display all columns

pool = index_components('000050.XSHG')
window = 3  # use the past 6 months for training and testing
risk_preference = -0.2  # in the range of (-1,1) exclusive!
verbose = True  # set to False when you don't want to waste time on plotting

today = datetime.today()
s = pool[11]
price_df = get_price(s,today-relativedelta(months=window+6),today-relativedelta(days=1),'1d',['open','close','high','low','volume'])
price_df['mkt'] = get_price('000300.XSHG',today-relativedelta(months=window+6),today-relativedelta(days=1),'1d',['close'])
X_ = pd.DataFrame(technical(price_df),index=price_df.index,columns=np.arange(78).astype('str')).ix[-window*20-1:,:]

for f in range(17,78): X_.ix[:,f] = X_.ix[:,f].astype('category')
y_ = (price_df.ix[X_.index,'close'] > price_df.ix[X_.index,'open']).astype('int')

X = X_.ix[:-1,:]
y = y_[1:]
print(X.shape)
print(y.shape)
print(y.value_counts())

y.describe()

X.describe(include=['number'])

X.describe(include=['category'])

unbalance = sum(y==1)/len(y)
if not (2/5 < unbalance < 3/5): print('Unbalanced dataset!!!')
unbalance

data = X
data.loc[:,'y'] = y.values
data.columns

if verbose:
    ax = plt.axes()
    sns.countplot(x='y', data=data, ax=ax)
    ax.set_title('Actual y Distribution (total: {})'.format(len(data)))
    plt.show()

if verbose:
    f, axarr = plt.subplots(4, 4, figsize=(15, 15))

    for i in range(16):
        sns.boxplot(x=str(i), y='y', data=data, showmeans=True, ax=axarr[i//4,i%4])
        axarr[i//4,i%4].set_title('factor '+str(i))
    
    plt.tight_layout()
    plt.show()

if verbose:
    sns.pairplot(data, vars=[str(i) for i in range(1,17)], hue='y', size=3)
    plt.show()

for i in range(17,78):
    dummies = pd.get_dummies(data.ix[:,i])
    dummies = dummies.add_prefix("{}#".format(i))
    data.drop(str(i), axis=1, inplace=True)
    data = data.join(dummies)

mask = data.std() == 0
data_valid = data.ix[:,~mask]

X = data_valid.drop('y', axis=1)
y = data_valid['y']
print(X.columns.tolist())

vt = VarianceThreshold().fit(X)
ranks = (-vt.variances_).argsort().argsort()
feat_var_threshold = X.columns[ranks<=20]
print(feat_var_threshold.tolist())

model = RandomForestClassifier()
model.fit(X,y)
feature_imp = pd.DataFrame(model.feature_importances_, index=X.columns, columns=['importance'])
feat_imp_20 = feature_imp.sort_values('importance', ascending=False).head(20).index
print(feat_imp_20.tolist())

X_minmax = MinMaxScaler([0,1]).fit_transform(X)
X_scored = SelectKBest(score_func=chi2, k='all').fit(X_minmax, y)
feature_scoring = pd.DataFrame({'feature': X.columns,'score': X_scored.scores_})
feat_scored_20 = feature_scoring.sort_values('score', ascending=False).head(20)['feature'].values
print(feat_scored_20.tolist())

rfe = RFE(LogisticRegression(),20)
rfe.fit(X, y)
feature_rfe_scoring = pd.DataFrame({'feature': X.columns,'score': rfe.ranking_})
feat_rfe_20 = feature_rfe_scoring[feature_rfe_scoring['score'] == 1]['feature'].values
print(feat_rfe_20.tolist())

features = np.hstack([feat_var_threshold,feat_imp_20,feat_scored_20,feat_rfe_20])
features = np.unique(features)
print('Final features ({} in total): {}'.format(len(features),', '.join(features)))

data_clean = data_valid.ix[:,features.tolist()+['y']]
tot_len = len(data_clean)
split_len = tot_len * 3 // 4
X_train = data_clean.ix[:split_len,features]
X_test = data_clean.ix[split_len:,features]
y_train = data_clean.ix[:split_len,'y']
y_test = data_clean.ix[split_len:,'y']
print('Clean dataset shape: {}'.format(data_clean.shape))
print('Train features shape: {}'.format(X_train.shape))
print('Test features shape: {}'.format(X_test.shape))
print('Train label shape: {}'. format(y_train.shape))
print('Test label shape: {}'. format(y_test.shape))

if verbose:
    components = 8
    pca = PCA(n_components=components).fit(X_train)
    pca_variance_explained_df = pd.DataFrame({'component': np.arange(1, components+1),'variance_explained': pca.explained_variance_ratio_})
    ax = sns.barplot(x='component', y='variance_explained', data=pca_variance_explained_df)
    ax.set_title('PCA - Variance explained')
    plt.show()

if verbose:
    X_pca = pd.DataFrame(pca.transform(X_train)[:,:2])
    X_pca['target'] = y_train.values
    X_pca.columns = ["x", "y", "target"]
    sns.lmplot('x','y', 
               data=X_pca, 
               hue="target", 
               fit_reg=False, 
               markers=["o", "x"], 
               palette="Set1", 
               size=7,
               scatter_kws={"alpha": .2})
    plt.show()

f = lambda x, r: (2*(x>0.5)-1)*abs(2*x-1)**(2**logit((r+1)/2))

if verbose:
    colors = ['red','orange','yellow','lime','cyan']
    r_grid = [-0.9,-0.5,0,0.5,0.9]
    x_grid = np.linspace(0,1,100)
    fig, ax = plt.subplots(figsize=(5,5))
    for i in range(len(r_grid)):
        f_grid = [f(x,r_grid[i]) for x in x_grid]
        ax.plot(x_grid, f_grid,c=colors[i],label='r = {:> 3.1f}'.format(r_grid[i]))
    ax.plot((0,1),(0,0),c='black')
    box = ax.get_position()
    ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
    ax.set_ylim([-1,1])
    ax.set_xlim([0,1])
    plt.xlabel('Accuracy')
    plt.ylabel('Utility')
    plt.title('Utility Curve for Different Risk Preferences')
    plt.show()

def custom_score(y_true, y_pred):
    x = sum(y_pred == y_true) / len(y_true)
    r = risk_preference
    return f(x,r)

seed = 7
processors = 1
num_folds = 4
num_instances = len(X_train)
scoring = make_scorer(custom_score)

kfold = KFold(n=num_instances, n_folds=num_folds, random_state=seed)

# Some basic models
models = []
models.append(('LR', LogisticRegression()))
models.append(('LDA', LinearDiscriminantAnalysis()))
models.append(('KNN', KNeighborsClassifier(n_neighbors=5)))
models.append(('DT', DecisionTreeClassifier()))
models.append(('GNB', GaussianNB()))
models.append(('SVC', SVC(probability=True, class_weight='balanced')))

# Evaluate each model in turn
results = []
names = []

for name, model in models:
    cv_results = cross_val_score(model, X_train, y_train, cv=kfold, scoring=scoring, n_jobs=processors)
    results.append(cv_results)
    names.append(name)
    print("{}: ({:.3f}) +/- ({:.3f})".format(name, cv_results.mean(), cv_results.std()))

cart = DecisionTreeClassifier()
num_trees = 100
model = BaggingClassifier(base_estimator=cart, n_estimators=num_trees, random_state=seed)
results = cross_val_score(model, X_train, y_train, cv=kfold, scoring=scoring, n_jobs=processors)
print("({:.3f}) +/- ({:.3f})".format(results.mean(), results.std()))

num_trees = 100
num_features = 10
model = RandomForestClassifier(n_estimators=num_trees, max_features=num_features)
results = cross_val_score(model, X_train, y_train, cv=kfold, scoring=scoring, n_jobs=processors)
print("({:.3f}) +/- ({:.3f})".format(results.mean(), results.std()))

num_trees = 100
num_features = 10
model = ExtraTreesClassifier(n_estimators=num_trees, max_features=num_features)
results = cross_val_score(model, X_train, y_train, cv=kfold, scoring=scoring, n_jobs=processors)
print("({:.3f}) +/- ({:.3f})".format(results.mean(), results.std()))

model = AdaBoostClassifier(n_estimators=100, random_state=seed)
results = cross_val_score(model, X_train, y_train, cv=kfold, scoring=scoring, n_jobs=processors)
print("({:.3f}) +/- ({:.3f})".format(results.mean(), results.std()))

model = GradientBoostingClassifier(n_estimators=100, random_state=seed)
results = cross_val_score(model, X_train, y_train, cv=kfold, scoring=scoring, n_jobs=processors)
print("({:.3f}) +/- ({:.3f})".format(results.mean(), results.std()))

model = XGBClassifier(n_estimators=100, seed=seed)
results = cross_val_score(model, X_train, y_train, cv=kfold, scoring=scoring, n_jobs=processors)
print("({:.3f}) +/- ({:.3f})".format(results.mean(), results.std()))

estimator_list = []
estimator_name_list = []
best_score_list = []
best_params_list = []

lr_grid = GridSearchCV(estimator = LogisticRegression(random_state=seed),
                       param_grid = {'penalty': ['l1', 'l2'],
                                     'C': [1e-3, 1e-2, 1, 10, 100, 1000]},
                       cv = kfold, 
                       scoring = scoring, 
                       n_jobs = processors)
lr_grid.fit(X_train, y_train)
print(lr_grid.best_score_)
print(lr_grid.best_params_)
estimator_name_list.append('LR')
estimator_list.append(LogisticRegression)
best_score_list.append(lr_grid.best_score_)
best_params_list.append(lr_grid.best_params_)

lda_grid = GridSearchCV(estimator = LinearDiscriminantAnalysis(),
                        param_grid = {'solver': ['svd','lsqr'],
                                      'n_components': [None, 2, 5, 10, 20]}, 
                        cv = kfold,
                        scoring = scoring, 
                        n_jobs = processors)
lda_grid.fit(X_train, y_train)
print(lda_grid.best_score_)
print(lda_grid.best_params_)
estimator_name_list.append('LDA')
estimator_list.append(LinearDiscriminantAnalysis)
best_score_list.append(lda_grid.best_score_)
best_params_list.append(lda_grid.best_params_)

dt_grid = GridSearchCV(estimator = DecisionTreeClassifier(random_state=seed),
                       param_grid = {'criterion': ['gini','entropy'],
                                     'max_depth': [None, 2, 5, 10], 
                                     'max_features': [None,'auto',20]},
                       cv = kfold,
                       scoring = scoring, 
                       n_jobs = processors)
dt_grid.fit(X_train, y_train)
print(dt_grid.best_score_)
print(dt_grid.best_params_)
estimator_name_list.append('DT')
estimator_list.append(DecisionTreeClassifier)
best_score_list.append(dt_grid.best_score_)
best_params_list.append(dt_grid.best_params_)

knn_grid = GridSearchCV(estimator = KNeighborsClassifier(),
                        param_grid = {'n_neighbors': [5, 10, 25],
                                      'algorithm': ['ball_tree'],
                                      'leaf_size': [2, 3, 5],
                                      'p': [1]}, 
                        cv = kfold, 
                        scoring = scoring, 
                        n_jobs = processors)
knn_grid.fit(X_train, y_train)
print(knn_grid.best_score_)
print(knn_grid.best_params_)
estimator_name_list.append('KNN')
estimator_list.append(KNeighborsClassifier)
best_score_list.append(knn_grid.best_score_)
best_params_list.append(knn_grid.best_params_)

rf_grid = GridSearchCV(estimator = RandomForestClassifier(warm_start=True, random_state=seed),
                       param_grid = {'n_estimators': [100, 200],
                                     'criterion': ['gini','entropy'],
                                     'max_features': [None, 20],
                                     'max_depth': [3, 5, 10],
                                     'bootstrap': [True]}, 
                       cv = kfold, 
                       scoring = scoring, 
                       n_jobs = processors)
rf_grid.fit(X_train, y_train)
print(rf_grid.best_score_)
print(rf_grid.best_params_)
estimator_name_list.append('RF')
estimator_list.append(RandomForestClassifier)
best_score_list.append(rf_grid.best_score_)
best_params_list.append(rf_grid.best_params_)

ext_grid = GridSearchCV(estimator = ExtraTreesClassifier(warm_start=True, random_state=seed),
                        param_grid = {'n_estimators': [100, 200],
                                      'criterion': ['gini', 'entropy'],
                                      'max_features': [None, 20],
                                      'max_depth': [3, 5, 10],
                                      'bootstrap': [True]}, 
                        cv = kfold, 
                        scoring = scoring, 
                        n_jobs = processors)
ext_grid.fit(X_train, y_train)
print(ext_grid.best_score_)
print(ext_grid.best_params_)
estimator_name_list.append('EXT')
estimator_list.append(ExtraTreesClassifier)
best_score_list.append(ext_grid.best_score_)
best_params_list.append(ext_grid.best_params_)

ada_grid = GridSearchCV(estimator = AdaBoostClassifier(random_state=seed),
                        param_grid = {'algorithm': ['SAMME', 'SAMME.R'],
                                      'n_estimators': [100, 200],
                                      'learning_rate': [1e-2, 1e-1]}, 
                        cv = kfold, 
                        scoring = scoring, 
                        n_jobs = processors)
ada_grid.fit(X_train, y_train)
print(ada_grid.best_score_)
print(ada_grid.best_params_)
estimator_name_list.append('ADA')
estimator_list.append(AdaBoostClassifier)
best_score_list.append(ada_grid.best_score_)
best_params_list.append(ada_grid.best_params_)

gbm_grid = GridSearchCV(estimator = GradientBoostingClassifier(warm_start=True, random_state=seed),
                        param_grid = {'n_estimators': [100, 200],
                                      'max_depth': [3, 5, 10],
                                      'max_features': [None, 20],
                                      'learning_rate': [1e-2, 1e-1]}, 
                        cv = kfold, 
                        scoring = scoring, 
                        n_jobs = processors)
gbm_grid.fit(X_train, y_train)
print(gbm_grid.best_score_)
print(gbm_grid.best_params_)
estimator_name_list.append('GBM')
estimator_list.append(GradientBoostingClassifier)
best_score_list.append(gbm_grid.best_score_)
best_params_list.append(gbm_grid.best_params_)

xgb_grid = GridSearchCV(estimator = XGBClassifier(nthread=1, seed=seed),
                        param_grid = {'n_estimators': [100, 200],
                                      'max_depth': [3, 5, 10],
                                      'min_child_weight': [1, 5],
                                      'gamma': [0, 0.1, 1, 10],
                                      'learning_rate': [1e-2, 1e-1]}, 
                        cv = kfold, 
                        scoring = scoring, 
                        n_jobs = processors)
xgb_grid.fit(X_train, y_train)
print(xgb_grid.best_score_)
print(xgb_grid.best_params_)
estimator_name_list.append('XGB')
estimator_list.append(XGBClassifier)
best_score_list.append(xgb_grid.best_score_)
best_params_list.append(xgb_grid.best_params_)

svc_grid = GridSearchCV(estimator = SVC(probability=True, class_weight='balanced'),
                        param_grid = {'C': [1e-3, 1e-2, 1e-1, 1, 10],
                                      'gamma': [1e-2, 1e-1, 1]},
                        cv = kfold, 
                        scoring = scoring, 
                        n_jobs = processors)
svc_grid.fit(X_train, y_train)
print(svc_grid.best_score_)
print(svc_grid.best_params_)
estimator_name_list.append('SVC')
estimator_list.append(SVC)
best_score_list.append(svc_grid.best_score_)
best_params_list.append(svc_grid.best_params_)

best_score_list_rounded = [round(s,3) for s in best_score_list]
seq = sorted(best_score_list, reverse=True)
rank = [seq.index(s) for s in best_score_list]

res = pd.DataFrame({'model':estimator_name_list,
                    'score':best_score_list_rounded,
                    'rank':rank},columns=['model','score','rank']).set_index('model')
res['rank'] = res['rank'].astype('object')
res.T

# Create sub models
estimators = []

# Choose models with larger logloss scores
for i in range(10):  # in total 10 models
    if rank[i] < 5:  # if in top 5
        estimators.append((estimator_name_list[i], estimator_list[i](**best_params_list[i])))

# Create the ensemble model
ensemble = VotingClassifier(estimators, voting='hard')

results = cross_val_score(ensemble, X_train, y_train, cv=kfold, scoring=scoring,n_jobs=processors)
print("({:.3f}) +/- ({:.3f})".format(results.mean(), results.std()))

model = ensemble
model.fit(X_train, y_train)
y_pred = model.predict(X_test)

print('Unbalance of the data: {:.3f}'.format(unbalance))

ac = accuracy_score(y_test, y_pred)
print('Accuracy: {:.3f}'.format(ac))  # (TP+TN) / (TP+FP+TN+FN)

pc = precision_score(y_test, y_pred)
print('Precision: {:.3f}'.format(pc))  # TP / (TP+FP)

re = recall_score(y_test, y_pred)
print('Recall: {:.3f}'.format(re))  # TP / (TP+FN)

cs = custom_score(y_test, y_pred)
print('Utility: {:.3f}'.format(cs))