{"cells":[{"metadata":{"_uuid":"5b12971a4ea735d69b63ef91428d867ac4968fb6"},"cell_type":"markdown","source":"\n**Hi, how are you?**\nI'm a big advocate of using algorithms from other fields (in this case finnace) for uncanny applications (in this case to predict earthquakes). \nI would like to put a big shout out to this git hub : https://github.com/borisbanushev/stockpredictionai\nWhich is amazing and from which I took out most of the code to adapt it to this problem, and you should definitely check that repo out !"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom utils import *\n\nimport os\nimport time\nimport numpy as np\n\nfrom mxnet import nd, autograd, gluon\nfrom mxnet.gluon import nn, rnn\nimport mxnet as mx\nimport datetime\nimport seaborn as sns\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom sklearn.decomposition import PCA\n\nimport math\n\nfrom sklearn.preprocessing import MinMaxScaler\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.preprocessing import StandardScaler\n\nimport xgboost as xgb\nfrom sklearn.metrics import accuracy_score\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\nfrom tqdm import tqdm_notebook as tqdm\nimport os\nimport gc\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9039e7023fcfcbef30d2a4c80031fbf90218a18a"},"cell_type":"markdown","source":"The input data"},{"metadata":{"trusted":true,"_uuid":"ee34e59813b05323f0a975a3fb2fd8639b29fb2f"},"cell_type":"code","source":"#This method is widely used by other kernels to implement data because the data is too large\n\ntrain = pd.read_csv('../input/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bf89b3919008152d7e2a491a11c1dad9f7d93b54"},"cell_type":"code","source":"#This method is widely used by other kernels to seperate the data.\n\nrows = 150_000\nsegments = int(np.floor(train.shape[0] / rows))\n\nx_train = pd.DataFrame(index=range(segments), dtype=np.int16,\n                       columns=['acoustic_data'])\n                       \ny_train = pd.DataFrame(index=range(segments), dtype=np.float32, columns=['time_to_failure'])\n                       \nfor segment in tqdm(range(segments)):\n    seg = train.iloc[segment*rows:segment*rows+rows]\n    x = seg['acoustic_data'].values\n    y = seg['time_to_failure'].values[-1]\n    y_train.loc[segment, 'time_to_failure'] = y\n    x_train.loc[segment, 'acoustic_data']=x.mean()                   ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bfe653f2a47c6fabfce926b88fa47fa81c16fe8e"},"cell_type":"markdown","source":"**Make some space**"},{"metadata":{"trusted":true,"_uuid":"a5e45925551a2d0c550c33b8925e0c50c34825da"},"cell_type":"code","source":"#My god this is 10 GB guys !! What the hell, delete it\ndel train\n#Collect garbage\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d6bedd94d2db5f3ca03f28058a46250b581d6916"},"cell_type":"markdown","source":"Small plot to see what it looks like"},{"metadata":{"trusted":true,"_uuid":"192f18beeb82549f391c2063e8d982a8ee6b10da"},"cell_type":"code","source":"plt.plot(x_train['acoustic_data'][0:35])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a0ac8f0ad07bf151a5feb6c37ddec02306eee8a1"},"cell_type":"markdown","source":"### Technical Analysis\n**Technical indicators** are a big part of financial analysis, this computes them."},{"metadata":{"trusted":true,"_uuid":"a298bd33268ee3b9c94de63a064273b832e88f5e"},"cell_type":"code","source":"def get_technical_indicators(dataset):\n    # Create 7 and 21 days Moving Average\n    dataset['ma7'] = dataset['acoustic_data'].rolling(window=7).mean()\n    dataset['ma21'] = dataset['acoustic_data'].rolling(window=21).mean()\n    \n    # Create MACD\n    dataset['26ema'] = dataset['acoustic_data'].ewm(span=26).mean() #pd.ewma(dataset['Close'], span=26)\n    dataset['12ema'] = dataset['acoustic_data'].ewm(span=12).mean() #pd.ewma(dataset['Close'], span=12)\n    dataset['MACD'] = (dataset['12ema']-dataset['26ema'])\n\n    # Create Bollinger Bands\n    dataset['20sd'] = dataset['acoustic_data'].rolling(20).std() #pd.stats.moments.rolling_std(dataset['Close'],20)\n    dataset['upper_band'] = dataset['ma21'] + (dataset['20sd']*2)\n    dataset['lower_band'] = dataset['ma21'] - (dataset['20sd']*2)\n    \n    # Create Exponential moving average\n    dataset['ema'] = dataset['acoustic_data'].ewm(com=0.5).mean()\n    \n    # Create Momentum\n    dataset['momentum'] = dataset['acoustic_data']-1\n    \n    return dataset","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9a6f4db8583df7664703d9ab080192a0aacc3d7f"},"cell_type":"code","source":"\n\n#Lets get these technical indicators\ntechnical_train = get_technical_indicators(x_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c5f4c640283232c8a749eb6a9079bd00456e7f1f"},"cell_type":"code","source":"#Drop the NA values got from the rolling averages\ntechnical_train_na = technical_train.dropna()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3adf3da57b59555c11f3adadffd9c197441e5130"},"cell_type":"code","source":"def plot_technical_indicators(dataset, last_days):\n    plt.figure(figsize=(16, 10), dpi=100)\n    shape_0 = dataset.shape[0]\n    xmacd_ = shape_0-last_days\n    \n    #dataset = dataset.iloc[-last_days:, :]\n    x_ = range(3, dataset.shape[0])\n    x_ =list(dataset.index)\n    \n    # Plot first subplot\n    #plt.subplot(2, 1, 1)\n    plt.plot(dataset['ma7'],label='MA 7', color='g',linestyle='--')\n    plt.plot(dataset['acoustic_data'],label='Real Data', color='b')\n    plt.plot(dataset['ma21'],label='MA 21', color='r',linestyle='--')\n    plt.plot(dataset['upper_band'],label='Upper Band', color='c')\n    plt.plot(dataset['lower_band'],label='Lower Band', color='c')\n    plt.plot((y_train['time_to_failure']/10)+4,label='Time to failure',color = 'y')\n    plt.fill_between(x_, dataset['lower_band'], dataset['upper_band'], alpha=0.35)\n    plt.title('Technical indicators')\n    plt.ylabel('Signal')\n    plt.legend()\n\n    \n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"699c08bf954796607f0a27977b6b5c74b659d19b"},"cell_type":"markdown","source":"**Man does it look good ! **"},{"metadata":{"trusted":true,"_uuid":"92af643af1359f4acfb62121ebc9d14921315be0"},"cell_type":"code","source":"plot_technical_indicators(technical_train_na,1000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"029e638a042b810c9487647138ee1c1513d39151"},"cell_type":"code","source":"#Let's use the FFT to get more features.\ndata_FT = x_train[['acoustic_data']]\n\nclose_fft = np.fft.fft(np.asarray(data_FT['acoustic_data'].tolist()))\nfft_df = pd.DataFrame({'fft':close_fft})\nfft_df['absolute'] = fft_df['fft'].apply(lambda x: np.abs(x))\nfft_df['angle'] = fft_df['fft'].apply(lambda x: np.angle(x))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"79fc26b74f8ed9d15848c1c5fda0d9650ac453f6"},"cell_type":"code","source":"plt.figure(figsize=(14, 7), dpi=100)\nfft_list = np.asarray(fft_df['fft'].tolist())\nfor num_ in [3, 6, 9, 25, 100]:\n    fft_list_m10= np.copy(fft_list); fft_list_m10[num_:-num_]=0\n    plt.plot(np.fft.ifft(fft_list_m10), label='Fourier transform with {} components'.format(num_))\n    \n#plt.plot(data_FT['acoustic_data'],  label='Real')\nplt.xlabel('rows')\nplt.ylabel('amplitude')\nplt.title('Fourier transform')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"da9afa0f3f2b25391db89dd60149383bc6c99316"},"cell_type":"markdown","source":"**Here we use the Fourrier transform to get more features**"},{"metadata":{"trusted":true,"_uuid":"430109e36bdbf19e7924a641a13cf62b96e9cb1e"},"cell_type":"code","source":"from collections import deque\nitems = deque(np.asarray(fft_df['absolute'].tolist()))\nitems.rotate(int(np.floor(len(fft_df)/2)))\nplt.figure(figsize=(10, 7), dpi=80)\nplt.stem(items)\nplt.title('Components of Fourier transforms')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"57575ac152f702316f178680ad0a8ccd9c10c4f6"},"cell_type":"markdown","source":"**The ARIMA model** is big for stock prediction, let's see if it can also predict earthquake amplitude."},{"metadata":{"trusted":true,"_uuid":"fa3b5ce76c1c9e2308ec112aa9c0ca92547020a9"},"cell_type":"code","source":"from statsmodels.tsa.arima_model import ARIMA\nfrom pandas import DataFrame\nfrom pandas import datetime","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2408d74818dfa2bf0137790ea7f2b18a4746a8a5"},"cell_type":"code","source":"series = data_FT['acoustic_data']\nmodel = ARIMA(series, order=(5, 1, 0))\nmodel_fit = model.fit(disp=0)\nprint(model_fit.summary())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5de722f8c5664e881ed1de2ad3f5649424363c86"},"cell_type":"code","source":"from pandas.tools.plotting import autocorrelation_plot\nautocorrelation_plot(series)\nplt.figure(figsize=(10, 7), dpi=80)\nplt.show() ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8de2a5d6666b9972cc8e0630e33d30527cb7024c"},"cell_type":"code","source":"from pandas import read_csv\nfrom pandas import datetime\nfrom statsmodels.tsa.arima_model import ARIMA\nfrom sklearn.metrics import mean_squared_error\n\nX = series.values\nsize = int(len(X) * 0.01)\ntrain, test = X[0:size], X[size:len(X)]\nhistory = [x for x in train]\npredictions = list()\nfor t in range(len(test)):\n    model = ARIMA(history, order=(5,1,0))\n    model_fit = model.fit(disp=0)\n    output = model_fit.forecast()\n    yhat = output[0]\n    predictions.append(yhat)\n    obs = test[t]\n    history.append(obs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a015bb46f1b43e587f71355ad2e4ec2f50f44b68"},"cell_type":"code","source":"error = mean_squared_error(test, predictions)\nprint('Test MSE: %.3f' % error)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"43e155eadf0c7efae2f9e705cb2113413fbbc8e0"},"cell_type":"code","source":"plt.figure(figsize=(12, 6), dpi=100)\nplt.plot(test, label='Real')\nplt.plot(predictions, color='red', label='Predicted')\nplt.xlabel('Rows')\nplt.ylabel('Amplitude')\nplt.title('ARIMA model on Amplitude')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a9780bb92ec46ddac20e0627150e0cd2bdc0379f"},"cell_type":"code","source":"arima = [np.NaN for i in range(41)] + list(test) #complete list for total database concatenation\nARIMA = pd.DataFrame(data=arima, columns=['ARIMA'])\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f825dff831b08c788f3c0a2fe1ec393e81ac6a01"},"cell_type":"markdown","source":"**Wow looks good ! We'll try to use it in our model**"},{"metadata":{"trusted":true,"_uuid":"0650705b183858dcd89aa14a82969e4f972b970f"},"cell_type":"code","source":"df_total = pd.concat([technical_train,ARIMA,fft_df.drop(['fft'],axis=1),y_train],ignore_index=False,sort=False,axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6ff38aed968f917f954e736734d065b7c64234a2"},"cell_type":"code","source":"df_total.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"455904f07a2cdbf4244c5f7996dec3649a86c42c"},"cell_type":"code","source":"df_total_withoutna=df_total.dropna()\ndf_total_withoutna.iloc[:,:14].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"32c0b8bbb22d5fa33822fecdbaec8d698d4cddd2"},"cell_type":"code","source":"def get_feature_importance_data(data_income):\n    data = data_income.copy()\n    y = data['time_to_failure']\n    X = data.iloc[:, :14]\n    \n    train_samples = int(X.shape[0] * 0.65)\n \n    X_train = X.iloc[:train_samples]\n    X_test = X.iloc[train_samples:]\n\n    y_train = y.iloc[:train_samples]\n    y_test = y.iloc[train_samples:]\n    \n    return (X_train, y_train), (X_test, y_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dc8654b27ae74922923b2206ad6fdfcb5e89aad3"},"cell_type":"code","source":"# Get training and test data\n(X_train_FI, y_train_FI), (X_test_FI, y_test_FI) = get_feature_importance_data(df_total_withoutna)\nregressor = xgb.XGBRegressor(gamma=0.0,n_estimators=150,base_score=0.7,colsample_bytree=1,learning_rate=0.05)\nxgbModel = regressor.fit(X_train_FI,y_train_FI, \\\n                         eval_set = [(X_train_FI, y_train_FI), (X_test_FI, y_test_FI)], \\\n                         verbose=False)\neval_result = regressor.evals_result()\ntraining_rounds = range(len(eval_result['validation_0']['rmse']))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"241c7aa306e703fa5d47cee84a07c4757dafaef7"},"cell_type":"code","source":"plt.scatter(x=training_rounds,y=eval_result['validation_0']['rmse'],label='Training Error')\nplt.scatter(x=training_rounds,y=eval_result['validation_1']['rmse'],label='Validation Error')\nplt.xlabel('Iterations')\nplt.ylabel('RMSE')\nplt.title('Training Vs Validation Error')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"41e80d713eb5a0475bf29f02579067d84978b9b9"},"cell_type":"code","source":"fig = plt.figure(figsize=(8,8))\nplt.xticks(rotation='vertical')\nplt.bar([i for i in range(len(xgbModel.feature_importances_))], xgbModel.feature_importances_.tolist(), tick_label=X_test_FI.columns)\nplt.title('Feature importance of indicators.')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6560dca11c7d5734e5883ace9ad01b73ffd548a1"},"cell_type":"markdown","source":"**We drop** the features which are not useful"},{"metadata":{"trusted":true,"_uuid":"efa932806f3188ffb5604e050ea71a5be00e1f72"},"cell_type":"code","source":"df_final = df_total_withoutna.drop(['acoustic_data','12ema','ema','momentum','ARIMA'],axis=1)\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f6e376e40c78aedfc978f86c5703b57cabee2349"},"cell_type":"markdown","source":"**Let's build the network now that we have some features**"},{"metadata":{"trusted":true,"_uuid":"2987ca3895aec0c6a5ea6719b475bff486acbcc0"},"cell_type":"code","source":"import keras\nfrom keras.layers import Dense, Flatten\nfrom keras.layers import Conv1D, MaxPooling1D\nfrom keras.models import Model\nfrom keras.layers import Input\nfrom tensorflow.python.keras import Sequential\nfrom keras.wrappers.scikit_learn import KerasRegressor\nfrom sklearn.model_selection import cross_val_score\nfrom sklearn.model_selection import KFold\nfrom sklearn.preprocessing import StandardScaler","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5b367fd97d81bd8d2c0788e0503a32c726395ff7"},"cell_type":"code","source":"def get_feature_importance_data(data_income):\n    data = data_income.copy()\n    y = data['time_to_failure']\n    X = data.iloc[:, :9]\n    \n    train_samples = int(X.shape[0] * 0.65)\n \n    X_train = X.iloc[:train_samples]\n    X_test = X.iloc[train_samples:]\n\n    y_train = y.iloc[:train_samples]\n    y_test = y.iloc[train_samples:]\n    \n    return (X_train, y_train), (X_test, y_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b572d5387f0e7bb8a4a852d0b7499a59927ad3b2"},"cell_type":"code","source":"(X_train_FI, y_train_FI), (X_test_FI, y_test_FI) = get_feature_importance_data(df_final)\n# Get training and test data\nregressor = xgb.XGBRegressor(gamma=0.0,n_estimators=150,base_score=0.7,colsample_bytree=1,learning_rate=0.05)\nxgbModel = regressor.fit(X_train_FI,y_train_FI, \\\n                         eval_set = [(X_train_FI, y_train_FI), (X_test_FI, y_test_FI)], \\\n                         verbose=False)\neval_result = regressor.evals_result()\ntraining_rounds = range(len(eval_result['validation_0']['rmse']))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"83adfc3a422b54a46b0cb164862d6238ed9cdabc"},"cell_type":"code","source":"lol = regressor.predict(X_test_FI)\nplt.plot(lol-5)\nplt.plot(y_test_FI.values)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1f9df887335c8e6a2ed63cc25e309ddea9da84e9"},"cell_type":"markdown","source":"That's what XGBoost can do"},{"metadata":{"trusted":true,"_uuid":"be0cc3bd098e82d29b84c52c5d3de873c4b71941"},"cell_type":"code","source":"\nXt=X_train_FI.values\nYt=y_train_FI.values\n\nXt.shape\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"db86539a4002ad1e3ab1c3147245c47c9b0a60cb"},"cell_type":"code","source":"model = Sequential()\nmodel.add(Dense(32,input_shape = (9,),activation = 'relu'))\nmodel.add(Dense(32,activation = 'relu'))\nmodel.add(Dense(32,activation = 'relu'))\nmodel.add(Dense(1))\nmodel.compile(loss = 'mae',optimizer = 'adam')\nmodel.fit(Xt, Yt, epochs=300, verbose=0)\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"653783372cdba207a1617ca3125a33e87daf9092"},"cell_type":"code","source":"y_pred = model.predict(X_test_FI.values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"34fc355220d0d91c2156b800baf795d7cd6ab38f"},"cell_type":"code","source":"plt.plot((y_pred-5)*2.5,color='g')\nplt.plot(y_test_FI.values,color='r')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a9a30360534eb8e1ea8aadeeba1fdb3acf992319"},"cell_type":"markdown","source":"**That's** what a simple NN can do, hopefully we can find something better because it just looks like noise..."},{"metadata":{"_uuid":"b80761c951e835ed0ef7b516870df53e5778b3ae"},"cell_type":"markdown","source":"**SO**:\nBasically, this doesn't give good results, but anyway, I hope it was still useful to you and give you ideas to solve this problem, I just finished this and will update as time goes on, thanks for the read, please please please coment if you have ideas on how to use this and maybe we'll be able to work something out. Much luck to you on this project !"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}