{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport os\nimport time\nimport datetime\nimport gc\nimport seaborn as sns\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nfrom scipy import stats\nfrom tqdm import tqdm\nfrom tqdm import tqdm_notebook\n#timeseries libraries\nfrom statsmodels import tsa\nfrom statsmodels.graphics import tsaplots\nfrom statsmodels.stats.diagnostic import acorr_ljungbox\nfrom statsmodels.api import qqplot\nfrom statsmodels.tsa.seasonal import seasonal_decompose\nfrom statsmodels.tsa.stattools import adfuller\n\nfrom statsmodels.tsa.arima_model import ARMA \nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\n\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import mean_squared_error,mean_absolute_error","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"scrolled":false},"cell_type":"code","source":"%%time\ntrain = pd.read_csv('../input/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"97cf73bf1a9bc3e059d79ebcd2b85ca9dd24b166"},"cell_type":"code","source":"#test = pd.read_csv('../input/test/.csv')\n#test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a70e3c941105df14b0f06716ca955f91ee7a3cb7","scrolled":false},"cell_type":"code","source":"train_acoustic_data_small = train['acoustic_data'].values[::50]\ntrain_time_to_failure_small = train['time_to_failure'].values[::50]\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title(\"Trends of acoustic_data and time_to_failure. 2% of data (sampled)\")\nplt.plot(train_acoustic_data_small, color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(train_time_to_failure_small, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))\nplt.grid(False)\n\ndel train_acoustic_data_small\ndel train_time_to_failure_small","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6d90003ea6c5ed92d021fdae1f49815a4fd3856f"},"cell_type":"code","source":"train_time_to_failure_small = train['time_to_failure'].values[::50]\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nax2 = ax1.twinx()\nplt.plot(train_time_to_failure_small, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))\nplt.grid(False)\ndel train_time_to_failure_small","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"507a193f3dd43c32f28d111b6ee2160a4c882b3f"},"cell_type":"markdown","source":"# Preparing and exploring data using average and std devation."},{"metadata":{"trusted":true,"_uuid":"9abffe6895205b8176fb541525b8be0d22cc8d77","scrolled":true},"cell_type":"code","source":"#let's prepare our data\nseg_length = 60000\ntotal_samples = int(np.floor((train.shape[0]) / seg_length))\n\n#we will be using a total of nine different features as given below for making our predictions\ncols = ['average', 'std'] #our features used for the prediction\nx_train = pd.DataFrame(index = range(total_samples), columns = cols, dtype = np.float64) #an empty dataframe holding our feature values\ny_train = pd.DataFrame(index = range(total_samples), columns = ['time_to_failure'], dtype = np.float64) #an empty dataframe holding our target labels","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cbaba0a21f3ecd47b6cc40a8c3ad0f136c8c43f6"},"cell_type":"code","source":"for value in tqdm(range(total_samples)):\n    sample = train.iloc[value*seg_length : value*seg_length + seg_length]\n    x = sample['acoustic_data'].values\n    y = sample['time_to_failure'].values[-1]\n    \n    y_train.loc[value, 'time_to_failure'] = y\n    \n    x_train.loc[value, 'average'] = x.mean()\n    x_train.loc[value, 'std'] = x.std()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"97a3978bbbaf2b371b0e6d3726416de48c40f770"},"cell_type":"code","source":"train_acoustic_data_small = x_train['average']\ntrain_time_to_failure_small = y_train['time_to_failure']\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title(\"Trends of acoustic_data and time_to_failure. 2% of data (sampled)\")\nplt.plot(train_acoustic_data_small, color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(train_time_to_failure_small, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))\nplt.grid(False)\n\ndel train_acoustic_data_small\ndel train_time_to_failure_small","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"72b4c6f429fae21d7799f6f035b8d1846863b17e"},"cell_type":"code","source":"x_train['seismic'] = x_train.average.diff()\ntrain_seismic = x_train['seismic']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8bd4c3feb04a9fe590d869dd6b379b6d1033321d"},"cell_type":"code","source":"nlags=13  #define number of lags to plot on ACF/PACF plots\nfig = plt.figure(figsize=(12,8))\nax1 = fig.add_subplot(211)\nfig = tsaplots.plot_acf(train_seismic.dropna(), lags=nlags, ax=ax1,title='Autocorrelation of Seismic Data')\nplt.xticks(list(range(0,nlags)))\n\nax2 = fig.add_subplot(212)\nfig = tsaplots.plot_pacf(train_seismic.dropna(), lags=nlags, ax=ax2,title='Partial Autocorrelation of Seismic Data')\nplt.xticks(list(range(0,nlags)))\nplt.xlabel('Lag-k')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"be07da8e461550d062d0238dc148b4e488e41259"},"cell_type":"code","source":"train_seismic.dropna(inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"183e4df3827d1d77615c2982079532084f94d694"},"cell_type":"code","source":"fig=plt.figure(figsize=(16,6))\n\nax1=plt.subplot(121)\nfig=sns.distplot(train_seismic,bins=200,ax=ax1)\nplt.title('Distribution of Seismic Data')\n\nax2=plt.subplot(122)\nfig=  qqplot(train_seismic,fit=True,line='45',ax=ax2)\nplt.title('Initial Q-Q plot for Seismic Data')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"99e7eeefdaf10a493ed0e81fd0ee2c158f09e879"},"cell_type":"code","source":"#Jarque-Bera test for normality\nval,p=stats.jarque_bera(train_seismic)\n\nprint('Jarque-Bera Test Results:\\nStatistics = {}\\np-value = {}'.format(val,p))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"debd361a83c410f3ab64ead8e054f325bc0734d8"},"cell_type":"code","source":"# Ljung-box test for auto-correlation\n_,p=acorr_ljungbox(train_seismic,lags=10)\nprint('Ljung-Box test p-values for 10-lags:\\n',p)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"94f9bd009a810ddba458109e9cf38aff9dfba255"},"cell_type":"code","source":"test_stat,pval,usedlag,_,CI=adfuller(train_seismic,regression='c',autolag='BIC')[:5]\nprint('Augmented Dickey Fuller test results:')\nprint('Test Statistics:',test_stat)\nprint('p-value:',pval)\nprint('Used Lags:',usedlag)\nprint('Critical Values:')\nfor key, value in CI.items():\n    print('\\t%s: %.3f' % (key, value))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c5f78977d17f68c67f10c87a7a1ae88d72526237"},"cell_type":"markdown","source":"Null Hypothesis (H0): If accepted, it suggests the time series has a unit root, meaning it is non-stationary. It has some time dependent structure.   \nAlternate Hypothesis (H1): The null hypothesis is rejected; it suggests the time series does not have a unit root, meaning it is stationary. It does not have time-dependent structure. We interpret this result using the p-value from the test.    \nA p-value below a threshold (such as 5% or 1%) suggests we reject the null hypothesis (stationary), otherwise a p-value above the threshold suggests we accept the null hypothesis (non-stationary).  \np-value > 0.05: Accept the null hypothesis (H0), the data has a unit root and is non-stationary.     \np-value <= 0.05: Reject the null hypothesis (H0), the data does not have a unit root and is stationary.    \n\n\nRunning the test statistic value of -26.753. As per ADF statistic we can see that our statistic value of -26.753 is less than the value of -3.432 at 1%.\nThe p-value is < 0.05.\nComparing the p-value and test statistic to the critical values, it looks like we would have to Reject the null hypothesis (H0),  the data does not have a unit root and is stationary."},{"metadata":{"trusted":true,"_uuid":"4810b8b16398e50dd1fcef8610a1232f1c218074"},"cell_type":"code","source":"x_train.to_csv(\"x_train_105.csv\")\ny_train.to_csv(\"y_train_105.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d10cb5bca806fe916de40f7c888276ee2e7eedc9"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6839b37cc564028ce9cd4c1fad625da2dd3d96d2"},"cell_type":"markdown","source":"# Taking difference of the dataset and the preparing the data."},{"metadata":{"trusted":true,"_uuid":"351e8d513ec911839ab43adeda852e07a46279e9"},"cell_type":"code","source":"#let's prepare our data\nseg_length = 50000\ntotal_samples = int(np.floor((train.shape[0]) / seg_length))\n\n#we will be using a total of nine different features as given below for making our predictions\ncols = ['average', 'std'] #our features used for the prediction\nx_train = pd.DataFrame(index = range(total_samples), columns = cols, dtype = np.float64) #an empty dataframe holding our feature values\ny_train = pd.DataFrame(index = range(total_samples), columns = ['time_to_failure'], dtype = np.float64) #an empty dataframe holding our target labels","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f199c5dac1d204e7eab1fb65f8a38779faf4d9a7"},"cell_type":"code","source":"for value in tqdm(range(total_samples)):\n    sample = train.iloc[value*seg_length : value*seg_length + seg_length]\n    x = sample['acoustic_data'].values\n    y = sample['time_to_failure'].values[-1]\n    \n    y_train.loc[value, 'time_to_failure'] = y\n    \n    x_train.loc[value, 'average'] = x.mean()\n    x_train.loc[value, 'std'] = x.std()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"62f8f6f2b2e4e3d07ad174033f489b2e534944cd"},"cell_type":"code","source":"train_acoustic_data_small = x_train['average']\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title(\"Trends of acoustic_data and time_to_failure. 2% of data (sampled)\")\nplt.plot(train_acoustic_data_small, color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nplt.grid(False)\n\ndel train_acoustic_data_small","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5827efbba08e0fe4d8b76021aec86d16f91cc3d1"},"cell_type":"code","source":"y_train.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5c3dec3a6a19ba4acd1e4fe85f9b93bb269b2c74"},"cell_type":"code","source":"nlags=60  #define number of lags to plot on ACF/PACF plots\nfig = plt.figure(figsize=(12,8))\nax1 = fig.add_subplot(211)\nfig = tsaplots.plot_acf(y_train, lags=nlags, ax=ax1,title='Autocorrelation of Seismic Data')\nplt.xticks(list(range(0,nlags)))\n\nax2 = fig.add_subplot(212)\nfig = tsaplots.plot_pacf(y_train, lags=nlags, ax=ax2,title='Partial Autocorrelation of Seismic Data')\nplt.xticks(list(range(0,nlags)))\nplt.xlabel('Lag-k')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9e84cc21b4714df7a668b6997277e6bbdd47d664"},"cell_type":"markdown","source":"# Building a Time Series Regression Model"},{"metadata":{"trusted":true,"_uuid":"6fb871b6031120fffef204c09087957b20129ec1"},"cell_type":"code","source":"train_main = train.values[::6000]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"0243f2d18194687820275ed53fd30901a9a92945"},"cell_type":"code","source":"train_main = pd.DataFrame(train_main, columns=['seismic','ttf'])\ntrain_main.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"2a1bd47084e0ed78dc3e3410f5efb4a931fea665"},"cell_type":"code","source":"train_acoustic_data_small = train_main['seismic']\ntrain_time_to_failure_small = train_main['ttf']\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title(\"Trends of acoustic_data and time_to_failure. 2% of data (sampled)\")\nplt.plot(train_acoustic_data_small, color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(train_time_to_failure_small, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))\nplt.grid(False)\n\ndel train_acoustic_data_small\ndel train_time_to_failure_small","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c2e87c5a61a5ccd62df5af38d2a1a394fdfa8c66"},"cell_type":"code","source":"train_lr = train_main[['seismic']]\ntarget_lr = train_main[['ttf']]\nprint(train_lr.shape)\nprint(target_lr.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6657d14fcce8b342ed9c9f22625a521ecb3ef4d0"},"cell_type":"code","source":"plt.plot(train_main.ttf.diff())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"35003f8a79151032c2630e176dc8b8bd772fef4a"},"cell_type":"code","source":"lr = LinearRegression()\nmodel = lr.fit(train_lr,target_lr)\npredict = model.predict(train_lr)\nresiduals = model.predict(train_lr) - target_lr","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8543a9b7a3f4d51ec08a5ab035dfe27c6a09aa8f"},"cell_type":"code","source":"mean_absolute_error(target_lr, predict)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ff9ba077936d1dc48385142ed89b64c1a095d567"},"cell_type":"code","source":"'''fig = plt.figure(figsize=(16,10))\nnlags = 15\nax1 = plt.subplot(2, 2, 1)\nfig = tsaplots.plot_acf(residuals, lags=nlags, ax=ax1,title='ACF for residuals')\nplt.xticks(list(range(0,nlags)))\n\nax2 = plt.subplot(2, 2, 2)\nfig = tsaplots.plot_pacf(residuals, lags=nlags, ax=ax2,title='PACF for residuals')\nplt.xticks(list(range(0,nlags)))\nplt.xlabel('Lag-k')\nplt.show() '''","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"ddadf468f5fdc39443fbf47f8225b7060326f99b"},"cell_type":"code","source":"model = SARIMAX(train_main.ttf, order=(1,0,0), exog=train_main[['seismic']])\nresult = model.fit()\nresult.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"19d8a17432a61493018fd6751d4449f468f1225b"},"cell_type":"code","source":"sample_submission = pd.read_csv(\"../input/sample_submission.csv\")\nsample_submission.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a9767af4aa2b1f73ca4adf2960d432f0f06822ca"},"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv', index_col='seg_id')\nX_test = pd.DataFrame(columns=train_main.columns, dtype=np.float64, index=submission.index)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8f46c567470fa4065e5d7e168336a1716df32b3e"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"12067100dc9e27e869fabdcc8bb4d1d81da49cf9"},"cell_type":"code","source":"test_result = []","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"d90a945cf25fe9c47425f5b999b730b97cdde5d2"},"cell_type":"code","source":"for seg_id in tqdm(X_test.index):\n    test = pd.read_csv('../input/test/' + seg_id + '.csv')\n    pred_uc = result.get_forecast(steps = len(test), exog=test[['acoustic_data']])\n    test_result.append(max(pred_uc.predicted_mean))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"4ce16d7e80dea171dca1eae1e275dfda8faad819"},"cell_type":"code","source":"test_result_submission = pd.DataFrame(test_result).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"022148c71febec8dac4a94efd2fcbfb85ca1107c"},"cell_type":"code","source":"test_result_submission = test_result_submission.T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6a6542843d03bc23975af296427778a3d68d633b"},"cell_type":"code","source":"test_result_submission.to_csv(\"submission\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5b52be4ff45be7c1a92def6a75c0c4dd7a23066e"},"cell_type":"code","source":"print(test_result_submission[:50])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7e11cea2273405f7b79681af0d0cf826d0722dfa"},"cell_type":"code","source":"print(test_result_submission[50:100])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6e7e06885742c7f2efb84160a0c70cb725286b89"},"cell_type":"code","source":"print(test_result_submission[100:150])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"641f0d7b11c6cc84d8c18ca875588472986905f5"},"cell_type":"code","source":"print(test_result_submission[150:200])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cf8a21455dd18d25871625ff16ef94f5f18fdac5"},"cell_type":"code","source":"print(test_result_submission[200:250])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"36e3198f9cf20c25ef24c1af057f1041f257fb32"},"cell_type":"code","source":"print(test_result_submission[250:300])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6c5369cda277c6b94bf8c296a02edc9bc58ed6a6"},"cell_type":"code","source":"print(test_result_submission[300:350])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f1666cff59de102456f9f70c2e084ebff6b932f1"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}