{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"<h2><center>Deep Learning for Time Series Forecasting</center></h2>\n\n<img src=\"https://raw.githubusercontent.com/dimitreOliveira/MachineLearning/master/Kaggle/Store%20Item%20Demand%20Forecasting%20Challenge/time-series%20graph.png\" width=\"800\">\n\n### The goal of this notebook is to develop and compare different approaches to time-series problems.","metadata":{"_uuid":"82b3b6c2dabd0ad00fcebe67ca5626753fb9365c"}},{"cell_type":"code","source":"import warnings\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom keras import optimizers\nfrom keras.utils import plot_model\nfrom keras.models import Sequential, Model\nfrom keras.layers.convolutional import Conv1D, MaxPooling1D\nfrom keras.layers import Dense, LSTM, SimpleRNN, RepeatVector, TimeDistributed, Flatten\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.model_selection import train_test_split\nimport plotly.plotly as py\nimport plotly.graph_objs as go\nfrom plotly.offline import init_notebook_mode, iplot\n\n# Basic packages\nimport datetime # manipulating date formats\nimport seaborn as sns # for prettier plots\n\n\n# TIME SERIES\nfrom statsmodels.tsa.arima_model import ARIMA\nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\nfrom pandas.plotting import autocorrelation_plot\nfrom statsmodels.tsa.stattools import adfuller, acf, pacf,arma_order_select_ic\nimport statsmodels.formula.api as smf\nimport statsmodels.tsa.api as smt\nimport statsmodels.api as sm\nimport scipy.stats as scs\n\n\n# settings\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\n%matplotlib inline\nwarnings.filterwarnings(\"ignore\")\ninit_notebook_mode(connected=True)\n\n# Set seeds to make the experiment more reproducible.\nfrom tensorflow import set_random_seed\nfrom numpy.random import seed\nset_random_seed(1)\nseed(1)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-08-12T15:08:37.224163Z","iopub.execute_input":"2022-08-12T15:08:37.224458Z","iopub.status.idle":"2022-08-12T15:08:37.284314Z","shell.execute_reply.started":"2022-08-12T15:08:37.224397Z","shell.execute_reply":"2022-08-12T15:08:37.283431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Loading data","metadata":{"_uuid":"dcc9f061fb8acd9447b60551e33f8f8d83fd9dad"}},{"cell_type":"code","source":"train = pd.read_csv('../input/demand-forecasting-kernels-only/train.csv', parse_dates=['date'])\n# read test data","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","execution":{"iopub.status.busy":"2022-08-12T14:42:35.111551Z","iopub.execute_input":"2022-08-12T14:42:35.111873Z","iopub.status.idle":"2022-08-12T14:42:35.935172Z","shell.execute_reply.started":"2022-08-12T14:42:35.111819Z","shell.execute_reply":"2022-08-12T14:42:35.934342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Train set","metadata":{"_uuid":"a60085dbb63091a2947a1f19bdbf2519b0370f1a"}},{"cell_type":"code","source":"# describe the dataframe ","metadata":{"_uuid":"ae977a0cb52dbbb8701bdf9c4d4393737a98d103","_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-12T12:01:09.311261Z","iopub.execute_input":"2022-08-12T12:01:09.311549Z","iopub.status.idle":"2022-08-12T12:01:09.360744Z","shell.execute_reply.started":"2022-08-12T12:01:09.311498Z","shell.execute_reply":"2022-08-12T12:01:09.360124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# show first 5 rows of the dataframe","metadata":{"_uuid":"9bb6d25dceea154555cabb6c45a83b94fb8dccf1","_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-12T11:25:25.765833Z","iopub.execute_input":"2022-08-12T11:25:25.766176Z","iopub.status.idle":"2022-08-12T11:25:25.781700Z","shell.execute_reply.started":"2022-08-12T11:25:25.766116Z","shell.execute_reply":"2022-08-12T11:25:25.780557Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Time period of the train dataset","metadata":{"_uuid":"827c9c89b0e144fd4cf4bf2e048e50f98cbe3e26"}},{"cell_type":"code","source":"# check the min, max date of the dataset","metadata":{"_uuid":"98d7855038ae46908b4dbdb9328ca500a6b5d552","_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-12T11:25:27.093761Z","iopub.execute_input":"2022-08-12T11:25:27.094076Z","iopub.status.idle":"2022-08-12T11:25:27.107778Z","shell.execute_reply.started":"2022-08-12T11:25:27.094000Z","shell.execute_reply":"2022-08-12T11:25:27.106957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Let's find out what's the time gap between the last day from training set from the last day of the test set, this will be out lag (the amount of day that need to be forecast)","metadata":{"_uuid":"f67cb93f552c3c7bac8cb8431072501348c4fa53"}},{"cell_type":"code","source":"print('Max date from train set: %s' % train['date'].max().date())\nprint('Max date from test set: %s' % test['date'].max().date())\n# calculate the lag_size between train and test dataset\n\nprint('Forecast lag size', lag_size)","metadata":{"_uuid":"7dd5441edc740bc3254694abb748f49b7ec3bbee","_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-12T11:25:28.926122Z","iopub.execute_input":"2022-08-12T11:25:28.926535Z","iopub.status.idle":"2022-08-12T11:25:28.946601Z","shell.execute_reply.started":"2022-08-12T11:25:28.926475Z","shell.execute_reply":"2022-08-12T11:25:28.945754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Basic EDA\n\nTo explore the time series data first we need to aggregate the sales by day","metadata":{"_uuid":"455b661280d3dff5b73f43df3fcd82e96e251663"}},{"cell_type":"code","source":"daily_sales = train.groupby('date', as_index=False)['sales'].sum()\nstore_daily_sales = train.groupby(['store', 'date'], as_index=False)['sales'].sum()\nitem_daily_sales = # calculate item daily sales","metadata":{"_uuid":"a865dab2049d6157f5c8b1bc90cc37c37b559854","_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-12T11:25:45.018830Z","iopub.execute_input":"2022-08-12T11:25:45.019199Z","iopub.status.idle":"2022-08-12T11:25:45.249653Z","shell.execute_reply.started":"2022-08-12T11:25:45.019129Z","shell.execute_reply":"2022-08-12T11:25:45.248835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Overall daily sales","metadata":{"_uuid":"153d5030e3381de2679042b1cd77a8889d616e47"}},{"cell_type":"code","source":"daily_sales_sc = go.Scatter(x=daily_sales[], y=daily_sales[])\nlayout = go.Layout(title='Daily sales', xaxis=dict(title='Date'), yaxis=dict(title='Sales'))\nfig = go.Figure(data=[daily_sales_sc], layout=layout)\niplot(fig)","metadata":{"_kg_hide-input":true,"_uuid":"79a4ccddbc65b201d6f7b4f55b152ef9eed6762f","execution":{"iopub.status.busy":"2022-08-12T11:25:46.304333Z","iopub.execute_input":"2022-08-12T11:25:46.304618Z","iopub.status.idle":"2022-08-12T11:25:47.498820Z","shell.execute_reply.started":"2022-08-12T11:25:46.304565Z","shell.execute_reply":"2022-08-12T11:25:47.498009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Daily sales by store","metadata":{"_uuid":"387e71579b0a4b5aca85a465c8e1b1a3a3289915"}},{"cell_type":"code","source":"store_daily_sales_sc = []\nfor store in store_daily_sales['store'].unique():\n    current_store_daily_sales = store_daily_sales[(store_daily_sales['store'] == store)]\n    store_daily_sales_sc.append(go.Scatter(x=current_store_daily_sales['date'], y=current_store_daily_sales['sales'], name=('Store %s' % store)))\n\nlayout = go.Layout(title='Store daily sales', xaxis=dict(title='Date'), yaxis=dict(title='Sales'))\nfig = go.Figure(data=store_daily_sales_sc, layout=layout)\niplot(fig)","metadata":{"_kg_hide-input":true,"_uuid":"0e8715f5557eac3cf3bd92eedec484c76ded1f61","execution":{"iopub.status.busy":"2022-08-12T11:25:48.197426Z","iopub.execute_input":"2022-08-12T11:25:48.197716Z","iopub.status.idle":"2022-08-12T11:25:49.106325Z","shell.execute_reply.started":"2022-08-12T11:25:48.197665Z","shell.execute_reply":"2022-08-12T11:25:49.105669Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Daily sales by item","metadata":{"_uuid":"15c99b3e7c8fdc1076ff4359697107c15d0da41a"}},{"cell_type":"code","source":"item_daily_sales_sc = []\nfor item in item_daily_sales['item'].unique():\n    current_item_daily_sales = item_daily_sales[(item_daily_sales['item'] == item)]\n    item_daily_sales_sc.append(go.Scatter(x=current_item_daily_sales['date'], y=current_item_daily_sales['sales'], name=('Item %s' % item)))\n\nlayout = go.Layout(title='Item daily sales', xaxis=dict(title='Date'), yaxis=dict(title='Sales'))\nfig = go.Figure(data=item_daily_sales_sc, layout=layout)\niplot(fig)","metadata":{"_kg_hide-input":true,"_uuid":"9036b41d4368d8ead697adbcafb50be91ca80c5c","execution":{"iopub.status.busy":"2022-08-12T11:26:00.457491Z","iopub.execute_input":"2022-08-12T11:26:00.457783Z","iopub.status.idle":"2022-08-12T11:26:05.019401Z","shell.execute_reply.started":"2022-08-12T11:26:00.457730Z","shell.execute_reply":"2022-08-12T11:26:05.018336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sales = pd.read_csv(\"../input/competitive-data-science-predict-future-sales/sales_train.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-08-12T14:48:00.144729Z","iopub.execute_input":"2022-08-12T14:48:00.145040Z","iopub.status.idle":"2022-08-12T14:48:02.850345Z","shell.execute_reply.started":"2022-08-12T14:48:00.144984Z","shell.execute_reply":"2022-08-12T14:48:02.849627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ts=sales.groupby([\"date_block_num\"])[\"item_cnt_day\"].sum()\nts.astype('float')\nplt.figure(figsize=(16,8))\nplt.title('Total Sales of the company')\nplt.xlabel('Time')\nplt.ylabel('Sales')\nplt.plot(ts);","metadata":{"execution":{"iopub.status.busy":"2022-08-12T14:48:41.030246Z","iopub.execute_input":"2022-08-12T14:48:41.030910Z","iopub.status.idle":"2022-08-12T14:48:41.598337Z","shell.execute_reply.started":"2022-08-12T14:48:41.030610Z","shell.execute_reply":"2022-08-12T14:48:41.597626Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(16,6))\nplt.plot(ts.rolling(window=12,center=False).mean(),label='Rolling Mean');\nplt.plot(ts.rolling(window=12,center=False).std(),label='Rolling sd');\nplt.legend();","metadata":{"execution":{"iopub.status.busy":"2022-08-12T15:02:46.500897Z","iopub.execute_input":"2022-08-12T15:02:46.501198Z","iopub.status.idle":"2022-08-12T15:02:46.895443Z","shell.execute_reply.started":"2022-08-12T15:02:46.501144Z","shell.execute_reply":"2022-08-12T15:02:46.894708Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import statsmodels.api as sm\n# multiplicative\nres = sm.tsa.seasonal_decompose(ts.values,freq=12,model=\"multiplicative\")\n#plt.figure(figsize=(16,12))\nfig = res.plot()\n#fig.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-12T15:06:20.409848Z","iopub.execute_input":"2022-08-12T15:06:20.410132Z","iopub.status.idle":"2022-08-12T15:06:20.993102Z","shell.execute_reply.started":"2022-08-12T15:06:20.410079Z","shell.execute_reply":"2022-08-12T15:06:20.992077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Additive model\nres = sm.tsa.seasonal_decompose(ts.values,freq=12,model=\"additive\")\n#plt.figure(figsize=(16,12))\nfig = res.plot()\n#fig.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-12T15:06:41.071076Z","iopub.execute_input":"2022-08-12T15:06:41.071369Z","iopub.status.idle":"2022-08-12T15:06:41.788160Z","shell.execute_reply.started":"2022-08-12T15:06:41.071313Z","shell.execute_reply":"2022-08-12T15:06:41.787010Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"we assume an additive model, then we can write\n\n> yt=St+Tt+Et \n\nwhere yt is the data at period t, St is the seasonal component at period t, Tt is the trend-cycle component at period tt and Et is the remainder (or irregular or error) component at period t\nSimilarly for Multiplicative model,\n\n> yt=St  x Tt x Et \n\n## Stationarity:\n\n![q](https://static1.squarespace.com/static/53ac905ee4b003339a856a1d/t/5818f84aebbd1ac01c275bac/1478031479192/?format=750w)\n\nStationarity refers to time-invariance of a series. (ie) Two points in a time series are related to each other by only how far apart they are, and not by the direction(forward/backward)\n\nWhen a time series is stationary, it can be easier to model. Statistical modeling methods assume or require the time series to be stationary.\n\n\nThere are multiple tests that can be used to check stationarity.\n* ADF( Augmented Dicky Fuller Test) \n* KPSS \n* PP (Phillips-Perron test)\n\nLet's just perform the ADF which is the most commonly used one.\n\nNote: [Step by step guide to perform dicky fuller test in Excel](http://www.real-statistics.com/time-series-analysis/stochastic-processes/dickey-fuller-test/)\n\n[Another Useful guide](http://www.blackarbs.com/blog/time-series-analysis-in-python-linear-models-to-garch/11/1/2016#AR) \n\n[good reference](https://github.com/ultimatist/ODSC17/blob/master/Time%20Series%20with%20Python%20(ODSC)%20STA.ipynb)\n","metadata":{}},{"cell_type":"code","source":"# Stationarity tests\ndef test_stationarity(timeseries):\n    \n    #Perform Dickey-Fuller test:\n    print('Results of Dickey-Fuller Test:')\n    dftest = adfuller(timeseries, autolag='AIC')\n    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])\n    for key,value in dftest[4].items():\n        dfoutput['Critical Value (%s)'%key] = value\n    print (dfoutput)\n\ntest_stationarity(ts)","metadata":{"execution":{"iopub.status.busy":"2022-08-12T15:08:40.635264Z","iopub.execute_input":"2022-08-12T15:08:40.635560Z","iopub.status.idle":"2022-08-12T15:08:40.659890Z","shell.execute_reply.started":"2022-08-12T15:08:40.635505Z","shell.execute_reply":"2022-08-12T15:08:40.659159Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# to remove trend\nfrom pandas import Series as Series\n# create a differenced series\ndef difference(dataset, interval=1):\n    diff = list()\n    for i in range(interval, len(dataset)):\n        value = dataset[i] - dataset[i - interval]\n        diff.append(value)\n    return Series(diff)\n\n# invert differenced forecast\ndef inverse_difference(last_ob, value):\n    return value + last_ob","metadata":{"execution":{"iopub.status.busy":"2022-08-12T15:09:11.416976Z","iopub.execute_input":"2022-08-12T15:09:11.417277Z","iopub.status.idle":"2022-08-12T15:09:11.422805Z","shell.execute_reply.started":"2022-08-12T15:09:11.417220Z","shell.execute_reply":"2022-08-12T15:09:11.421288Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ts=sales.groupby([\"date_block_num\"])[\"item_cnt_day\"].sum()\nts.astype('float')\nplt.figure(figsize=(16,16))\nplt.subplot(311)\nplt.title('Original')\nplt.xlabel('Time')\nplt.ylabel('Sales')\nplt.plot(ts)\nplt.subplot(312)\nplt.title('After De-trend')\nplt.xlabel('Time')\nplt.ylabel('Sales')\nnew_ts=difference(ts)\nplt.plot(new_ts)\nplt.plot()\n\nplt.subplot(313)\nplt.title('After De-seasonalization')\nplt.xlabel('Time')\nplt.ylabel('Sales')\nnew_ts=difference(ts,12)       # assuming the seasonality is 12 months long\nplt.plot(new_ts)\nplt.plot()","metadata":{"execution":{"iopub.status.busy":"2022-08-12T15:10:25.673866Z","iopub.execute_input":"2022-08-12T15:10:25.674158Z","iopub.status.idle":"2022-08-12T15:10:26.645404Z","shell.execute_reply.started":"2022-08-12T15:10:25.674105Z","shell.execute_reply":"2022-08-12T15:10:26.644560Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# now testing the stationarity again after de-seasonality\ntest_stationarity(new_ts)","metadata":{"execution":{"iopub.status.busy":"2022-08-12T15:14:35.963625Z","iopub.execute_input":"2022-08-12T15:14:35.963961Z","iopub.status.idle":"2022-08-12T15:14:35.980304Z","shell.execute_reply.started":"2022-08-12T15:14:35.963900Z","shell.execute_reply":"2022-08-12T15:14:35.979415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Sub-sample train set to get only the last year of data and reduce training time","metadata":{"_uuid":"dd74fe1b76da480ca147af0f68b1e07d479020fc"}},{"cell_type":"code","source":"train = train[(train['date'] >= '2017-01-01')]","metadata":{"_uuid":"58a105bc712d2e1da58dfa33cb62f1203de9f0a0","execution":{"iopub.status.busy":"2022-08-12T11:26:05.021267Z","iopub.execute_input":"2022-08-12T11:26:05.021790Z","iopub.status.idle":"2022-08-12T11:26:05.043044Z","shell.execute_reply.started":"2022-08-12T11:26:05.021726Z","shell.execute_reply":"2022-08-12T11:26:05.042416Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Rearrange dataset so we can apply shift methods","metadata":{"_uuid":"f30936a69f127157190ecfedcacfd3e50bbb1616"}},{"cell_type":"code","source":"train_gp = train.sort_values('date').groupby(['item', 'store', 'date'], as_index=False)\ntrain_gp = train_gp.agg({'sales':['mean']})\ntrain_gp.columns = ['item', 'store', 'date', 'sales']\ntrain_gp.head()","metadata":{"_uuid":"2f9ffc8f0946c5caf4d99dc878f83f9e29ee73dd","execution":{"iopub.status.busy":"2022-08-12T11:26:05.044133Z","iopub.execute_input":"2022-08-12T11:26:05.044549Z","iopub.status.idle":"2022-08-12T11:26:05.209560Z","shell.execute_reply.started":"2022-08-12T11:26:05.044505Z","shell.execute_reply":"2022-08-12T11:26:05.208810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Transform the data into a time series problem","metadata":{"_uuid":"72c037d182b0cd9b8d228d0042ace64c88e6cbd9"}},{"cell_type":"code","source":"def series_to_supervised(data, window=1, lag=1, dropnan=True):\n    cols, names = list(), list()\n    # Input sequence (t-n, ... t-1)\n    for i in range(window, 0, -1):\n        cols.append(data.shift(i))\n        names += [('%s(t-%d)' % (col, i)) for col in data.columns]\n    # Current timestep (t=0)\n    cols.append(data)\n    names += [('%s(t)' % (col)) for col in data.columns]\n    # Target timestep (t=lag)\n    cols.append(data.shift(-lag))\n    names += [('%s(t+%d)' % (col, lag)) for col in data.columns]\n    # Put it all together\n    agg = pd.concat(cols, axis=1)\n    agg.columns = names\n    # Drop rows with NaN values\n    if dropnan:\n        agg.dropna(inplace=True)\n    return agg","metadata":{"_uuid":"c9400d07223d75b50930e37194c2a682ca436bf2","execution":{"iopub.status.busy":"2022-08-12T11:26:05.210730Z","iopub.execute_input":"2022-08-12T11:26:05.211192Z","iopub.status.idle":"2022-08-12T11:26:05.218185Z","shell.execute_reply.started":"2022-08-12T11:26:05.211146Z","shell.execute_reply":"2022-08-12T11:26:05.217364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### We will use the current timestep and the last 29 to forecast 90 days ahead","metadata":{"_uuid":"76407684060139325ff834aa9215360c7f69b5a8"}},{"cell_type":"code","source":"window = 29\nlag = lag_size\nseries = series_to_supervised(train_gp.drop('date', axis=1), window=window, lag=lag)\nseries.head()","metadata":{"_uuid":"5cb604bd4fb5aa6a2fd966af6f7ed772a11ebf6e","execution":{"iopub.status.busy":"2022-08-12T11:26:05.511765Z","iopub.execute_input":"2022-08-12T11:26:05.512239Z","iopub.status.idle":"2022-08-12T11:26:06.374046Z","shell.execute_reply.started":"2022-08-12T11:26:05.512014Z","shell.execute_reply":"2022-08-12T11:26:06.373162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Drop rows with different item or store values than the shifted columns","metadata":{"_uuid":"f147b4d5ed1cde12436d07fb13fff6c074857f2f"}},{"cell_type":"code","source":"last_item = 'item(t-%d)' % window\nlast_store = 'store(t-%d)' % window\nseries = series[(series['store(t)'] == series[last_store])]\nseries = series[(series['item(t)'] == series[last_item])]","metadata":{"_uuid":"bcc650cb3df2f11a58b4288171663d67cab5f4f7","execution":{"iopub.status.busy":"2022-08-12T11:26:06.375511Z","iopub.execute_input":"2022-08-12T11:26:06.376027Z","iopub.status.idle":"2022-08-12T11:26:06.567596Z","shell.execute_reply.started":"2022-08-12T11:26:06.375811Z","shell.execute_reply":"2022-08-12T11:26:06.566873Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Remove unwanted columns","metadata":{"_uuid":"07c8a6ac6dbbe309b626d4e2f47a44fae987875c"}},{"cell_type":"code","source":"columns_to_drop = [('%s(t+%d)' % (col, lag)) for col in ['item', 'store']]\nfor i in range(window, 0, -1):\n    columns_to_drop += [('%s(t-%d)' % (col, i)) for col in ['item', 'store']]\nseries.drop(columns_to_drop, axis=1, inplace=True)\nseries.drop(['item(t)', 'store(t)'], axis=1, inplace=True)","metadata":{"_uuid":"830f58c2a435ce18271629008fb1cd79c9323ead","execution":{"iopub.status.busy":"2022-08-12T11:26:06.774168Z","iopub.execute_input":"2022-08-12T11:26:06.774426Z","iopub.status.idle":"2022-08-12T11:26:06.846041Z","shell.execute_reply.started":"2022-08-12T11:26:06.774377Z","shell.execute_reply":"2022-08-12T11:26:06.845243Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Train/validation split","metadata":{"_uuid":"cdb34762d10c594f403df21519abb22601bf4e01"}},{"cell_type":"code","source":"# Label\nlabels_col = 'sales(t+%d)' % lag_size\nlabels = series[labels_col]\nseries = series.drop(labels_col, axis=1)\n\nX_train, X_valid, Y_train, Y_valid = train_test_split(series, labels.values, test_size=0.4, random_state=0)\nprint('Train set shape', X_train.shape)\nprint('Validation set shape', X_valid.shape)\nX_train.head()","metadata":{"_uuid":"389bfd84bc41a330b2f5ec84dd0ca7517eb52e31","execution":{"iopub.status.busy":"2022-08-12T11:31:02.217927Z","iopub.execute_input":"2022-08-12T11:31:02.218257Z","iopub.status.idle":"2022-08-12T11:31:02.351684Z","shell.execute_reply.started":"2022-08-12T11:31:02.218196Z","shell.execute_reply":"2022-08-12T11:31:02.350820Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### MLP for Time Series Forecasting\n\n* First we will use a Multilayer Perceptron model or MLP model, here our model will have input features equal to the window size.\n* The thing with MLP models is that the model don't take the input as sequenced data, so for the model, it is just receiving inputs and don't treat them as sequenced data, that may be a problem since the model won't see the data with the sequence patter that it has.\n* Input shape **[samples, timesteps]**.","metadata":{"_uuid":"1c6dbf5e8699ea64e929009e7c47625eab207535"}},{"cell_type":"code","source":"epochs = \nbatch = \nlr = \nadam = optimizers.Adam(lr)","metadata":{"_uuid":"47962d535b4dddea0588a416ed742b23fdc44db0","execution":{"iopub.status.busy":"2022-08-12T11:31:02.611124Z","iopub.execute_input":"2022-08-12T11:31:02.611358Z","iopub.status.idle":"2022-08-12T11:31:02.633087Z","shell.execute_reply.started":"2022-08-12T11:31:02.611310Z","shell.execute_reply":"2022-08-12T11:31:02.632507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_mlp = Sequential()\nmodel_mlp.add(Dense(, activation='relu', input_dim=X_train.shape[1]))\nmodel_mlp.add(Dense())\nmodel_mlp.compile(loss=, optimizer=adam)\nmodel_mlp.summary()","metadata":{"_uuid":"da03328a583580969bd2603c1348b939d533092d","execution":{"iopub.status.busy":"2022-08-12T11:31:02.787414Z","iopub.execute_input":"2022-08-12T11:31:02.787632Z","iopub.status.idle":"2022-08-12T11:31:02.832062Z","shell.execute_reply.started":"2022-08-12T11:31:02.787588Z","shell.execute_reply":"2022-08-12T11:31:02.831200Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mlp_history = model_mlp.fit(X_train.values, Y_train, validation_data=(X_valid.values, Y_valid), epochs=epochs, verbose=2)","metadata":{"_uuid":"c741ddf5b382617173ef09450a836f940ab4f340","_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-08-12T11:31:03.472466Z","iopub.execute_input":"2022-08-12T11:31:03.472741Z","iopub.status.idle":"2022-08-12T11:36:25.573102Z","shell.execute_reply.started":"2022-08-12T11:31:03.472692Z","shell.execute_reply":"2022-08-12T11:36:25.572449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### CNN for Time Series Forecasting\n\n* For the CNN model we will use one convolutional hidden layer followed by a max pooling layer. The filter maps are then flattened before being interpreted by a Dense layer and outputting a prediction.\n* The convolutional layer should be able to identify patterns between the timesteps.\n* Input shape **[samples, timesteps, features]**.\n\n#### Data preprocess\n* Reshape from [samples, timesteps] into [samples, timesteps, features].\n* This same reshaped data will be used on the CNN and the LSTM model.","metadata":{"_uuid":"7ef4d50fedcf9516742b45f3e3a7dab2cb0d301b"}},{"cell_type":"code","source":"X_train_series = X_train.values.reshape((X_train.shape[0], X_train.shape[1], 1))\nX_valid_series = X_valid.values.reshape((X_valid.shape[0], X_valid.shape[1], 1))\nprint('Train set shape', X_train_series.shape)\nprint('Validation set shape', X_valid_series.shape)","metadata":{"_uuid":"00b13a9d49412dd9a50c7a213a042c47e05a53c0","execution":{"iopub.status.busy":"2022-08-12T11:45:57.754395Z","iopub.execute_input":"2022-08-12T11:45:57.754680Z","iopub.status.idle":"2022-08-12T11:45:57.789951Z","shell.execute_reply.started":"2022-08-12T11:45:57.754630Z","shell.execute_reply":"2022-08-12T11:45:57.788917Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_cnn = Sequential()\nmodel_cnn.add(Conv1D(filters=64, kernel_size=2, activation='relu', input_shape=(X_train_series.shape[1], X_train_series.shape[2])))\nmodel_cnn.add(MaxPooling1D(pool_size=2))\nmodel_cnn.add(Flatten())\nmodel_cnn.add(Dense(50, activation='relu'))\nmodel_cnn.add(Dense(1))\nmodel_cnn.compile(loss='mse', optimizer=adam)\nmodel_cnn.summary()","metadata":{"_uuid":"846fccfaa8d3761bbfc57b05ab447e9da07367ee","execution":{"iopub.status.busy":"2022-08-12T11:45:58.788828Z","iopub.execute_input":"2022-08-12T11:45:58.789361Z","iopub.status.idle":"2022-08-12T11:45:58.867889Z","shell.execute_reply.started":"2022-08-12T11:45:58.789302Z","shell.execute_reply":"2022-08-12T11:45:58.867290Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cnn_history = model_cnn.fit(X_train_series, Y_train, validation_data=(X_valid_series, Y_valid), epochs=epochs, verbose=2)","metadata":{"_uuid":"02870e84ad2997b21994732794251b5d5865211a","_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-08-12T11:46:00.131433Z","iopub.execute_input":"2022-08-12T11:46:00.131752Z","iopub.status.idle":"2022-08-12T11:53:03.410243Z","shell.execute_reply.started":"2022-08-12T11:46:00.131709Z","shell.execute_reply":"2022-08-12T11:53:03.409554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_rnn = Sequential()\nmodel_rnn.add(SimpleRNN(50, activation='relu', input_shape=(X_train_series.shape[1], X_train_series.shape[2])))\nmodel_rnn.add(Dense(1))\nmodel_rnn.compile(loss='mse', optimizer=adam)\nmodel_rnn.summary()","metadata":{"execution":{"iopub.status.busy":"2022-08-12T12:22:22.243268Z","iopub.execute_input":"2022-08-12T12:22:22.243562Z","iopub.status.idle":"2022-08-12T12:22:22.439271Z","shell.execute_reply.started":"2022-08-12T12:22:22.243509Z","shell.execute_reply":"2022-08-12T12:22:22.438604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### LSTM for Time Series Forecasting\n\n* Now the LSTM model actually sees the input data as a sequence, so it's able to learn patterns from sequenced data (assuming it exists) better than the other ones, especially patterns from long sequences.\n* Input shape **[samples, timesteps, features]**.","metadata":{"_uuid":"5eac02f27336cb62fa2de354ef647caa3fb4934e"}},{"cell_type":"code","source":"model_lstm = Sequential()\nmodel_lstm.add(LSTM(50, activation='relu', input_shape=(X_train_series.shape[1], X_train_series.shape[2])))\nmodel_lstm.add(Dense(1))\nmodel_lstm.compile(loss='mse', optimizer=adam)\nmodel_lstm.summary()","metadata":{"_uuid":"522965ca3e2490b8243032c6418ce50d5b724872","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"lstm_history = model_lstm.fit(X_train_series, Y_train, validation_data=(X_valid_series, Y_valid), epochs=15, verbose=2)","metadata":{"_uuid":"48ac362d2a96028dc2bba1477b3d357c4d79668d","_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### CNN-LSTM for Time Series Forecasting\n* Input shape **[samples, subsequences, timesteps, features]**.\n\n#### Model explanation from the [article](https://machinelearningmastery.com/how-to-get-started-with-deep-learning-for-time-series-forecasting-7-day-mini-course/)\n> \"The benefit of this model is that the model can support very long input sequences that can be read as blocks or subsequences by the CNN model, then pieced together by the LSTM model.\"\n>\n> \"When using a hybrid CNN-LSTM model, we will further divide each sample into further subsequences. The CNN model will interpret each sub-sequence and the LSTM will piece together the interpretations from the subsequences. As such, we will split each sample into 2 subsequences of 2 times per subsequence.\"\n>\n> \"The CNN will be defined to expect 2 timesteps per subsequence with one feature. The entire CNN model is then wrapped in TimeDistributed wrapper layers so that it can be applied to each subsequence in the sample. The results are then interpreted by the LSTM layer before the model outputs a prediction.\"\n\n#### Data preprocess\n* Reshape from [samples, timesteps, features] into [samples, subsequences, timesteps, features].","metadata":{"_uuid":"1ba1b9fede4966b07059698f25e7f18e58cab232"}},{"cell_type":"code","source":"subsequences = 2\ntimesteps = X_train_series.shape[1]//subsequences\nX_train_series_sub = X_train_series.reshape((X_train_series.shape[0], subsequences, timesteps, 1))\nX_valid_series_sub = X_valid_series.reshape((X_valid_series.shape[0], subsequences, timesteps, 1))\nprint('Train set shape', X_train_series_sub.shape)\nprint('Validation set shape', X_valid_series_sub.shape)","metadata":{"_uuid":"b971641bae57f404bf16dcb7fd4e6fcff47f2be7","execution":{"iopub.status.busy":"2022-08-12T12:22:35.038324Z","iopub.execute_input":"2022-08-12T12:22:35.038618Z","iopub.status.idle":"2022-08-12T12:22:35.045173Z","shell.execute_reply.started":"2022-08-12T12:22:35.038563Z","shell.execute_reply":"2022-08-12T12:22:35.044062Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_cnn_lstm = Sequential()\nmodel_cnn_lstm.add(TimeDistributed(Conv1D(filters=64, kernel_size=1, activation='relu'), input_shape=(None, X_train_series_sub.shape[2], X_train_series_sub.shape[3])))\nmodel_cnn_lstm.add(TimeDistributed(MaxPooling1D(pool_size=2)))\nmodel_cnn_lstm.add(TimeDistributed(Flatten()))\nmodel_cnn_lstm.add(LSTM(50, activation='relu'))\nmodel_cnn_lstm.add(Dense(1))\nmodel_cnn_lstm.compile(loss='mse', optimizer=adam)","metadata":{"_uuid":"57b15ea60310bc7503e76d553e2c6a0d68020c90","execution":{"iopub.status.busy":"2022-08-12T12:22:36.424095Z","iopub.execute_input":"2022-08-12T12:22:36.424388Z","iopub.status.idle":"2022-08-12T12:22:36.730766Z","shell.execute_reply.started":"2022-08-12T12:22:36.424337Z","shell.execute_reply":"2022-08-12T12:22:36.729964Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cnn_lstm_history = model_cnn_lstm.fit(X_train_series_sub, Y_train, validation_data=(X_valid_series_sub, Y_valid), epochs=epochs, verbose=2)","metadata":{"_uuid":"10c53d1c61050407ce9406093a0fb68d67895b27","_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Comparing models","metadata":{"_uuid":"624df546d7a4f95d8f9bd51d065e472568b25eb2","trusted":true}},{"cell_type":"code","source":"fig, axes = plt.subplots(2, 2, sharex=True, sharey=True,figsize=(22,12))\nax1, ax2 = axes[0]\nax3, ax4 = axes[1]\n\nax1.plot(mlp_history.history['loss'], label='Train loss')\nax1.plot(mlp_history.history['val_loss'], label='Validation loss')\nax1.legend(loc='best')\nax1.set_title('MLP')\nax1.set_xlabel('Epochs')\nax1.set_ylabel('MSE')\n\nax2.plot(cnn_history.history['loss'], label='Train loss')\nax2.plot(cnn_history.history['val_loss'], label='Validation loss')\nax2.legend(loc='best')\nax2.set_title('CNN')\nax2.set_xlabel('Epochs')\nax2.set_ylabel('MSE')\n\nax3.plot(lstm_history.history['loss'], label='Train loss')\nax3.plot(lstm_history.history['val_loss'], label='Validation loss')\nax3.legend(loc='best')\nax3.set_title('LSTM')\nax3.set_xlabel('Epochs')\nax3.set_ylabel('MSE')\n\nax4.plot(cnn_lstm_history.history['loss'], label='Train loss')\nax4.plot(cnn_lstm_history.history['val_loss'], label='Validation loss')\nax4.legend(loc='best')\nax4.set_title('CNN-LSTM')\nax4.set_xlabel('Epochs')\nax4.set_ylabel('MSE')\n\nplt.show()","metadata":{"_uuid":"47d797354baed8299917542f3a541e61e140e822","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### MLP on train and validation","metadata":{"_uuid":"13c7dc56802c04c698310f2fd424d2a001a3e995"}},{"cell_type":"code","source":"mlp_train_pred = model_mlp.predict(X_train.values)\nmlp_valid_pred = model_mlp.predict(X_valid.values)\nprint('Train rmse:', np.sqrt(mean_squared_error(Y_train, mlp_train_pred)))\nprint('Validation rmse:', np.sqrt(mean_squared_error(Y_valid, mlp_valid_pred)))","metadata":{"_uuid":"17866c9d0ae37cae8e0da41b1d93d7d4049ad043","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### CNN on train and validation","metadata":{"_uuid":"425f734d0e4c46059138987de98936899ec3bdae"}},{"cell_type":"code","source":"cnn_train_pred = model_cnn.predict(X_train_series)\ncnn_valid_pred = model_cnn.predict(X_valid_series)\nprint('Train rmse:', np.sqrt(mean_squared_error(Y_train, cnn_train_pred)))\nprint('Validation rmse:', np.sqrt(mean_squared_error(Y_valid, cnn_valid_pred)))","metadata":{"_uuid":"eaa3ca7385ee7c8d884f5e1e6abd54c0aea91293","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### LSTM on train and validation","metadata":{"_uuid":"2355313b77222db98a69564449942e3f8378336a"}},{"cell_type":"code","source":"lstm_train_pred = model_lstm.predict(X_train_series)\nlstm_valid_pred = model_cnn.predict(X_valid_series)\nprint('Train rmse:', np.sqrt(mean_squared_error(Y_train, lstm_train_pred)))\nprint('Validation rmse:', np.sqrt(mean_squared_error(Y_valid, lstm_valid_pred)))","metadata":{"_uuid":"939f86306a1325c949d8c97262147bcce4f88462","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### CNN-LSTM on train and validation","metadata":{"_uuid":"c48c38e1cf05cd7503126c5af9c460b743872b9d"}},{"cell_type":"code","source":"cnn_lstm_train_pred = model_cnn_lstm.predict(X_train_series_sub)\ncnn_lstm_valid_pred = model_cnn_lstm.predict(X_valid_series_sub)\nprint('Train rmse:', np.sqrt(mean_squared_error(Y_train, cnn_lstm_train_pred)))\nprint('Validation rmse:', np.sqrt(mean_squared_error(Y_valid, cnn_lstm_valid_pred)))","metadata":{"_uuid":"95e86d4884cb1c2769b194769f98fbc7da8987c3","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Conclusion\n\nHere you could see some approaches to a time series problem, how to develop and the differences between them, this is not meant to have a great performance, so if you want better results, you are more than welcomed to try a few different hyper-parameters, especially the window size and the networks topology, if you do, please let me know the results.\n\nI hope you learned a few things here, leave a feedback and if you liked what you saw make sure to check the [article](https://machinelearningmastery.com/how-to-get-started-with-deep-learning-for-time-series-forecasting-7-day-mini-course/) that I used as source.\n\nIf you want to check out how you can use LSTM as autoencoders and create new features that represent a time series take a look at my other kernel [Time-series forecasting with deep learning & LSTM autoencoders](https://www.kaggle.com/dimitreoliveira/time-series-forecasting-with-lstm-autoencoders/data).","metadata":{"_uuid":"adba93182e270328e63c5c0cb55d0d1ccb344660"}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}