{"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":"code","source":"import pandas as pd \nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport plotly.express as px\nimport plotly.graph_objects as go\nimport numpy as np\nfrom matplotlib.gridspec import GridSpec\nfrom fbprophet import Prophet\nfrom fbprophet.plot import plot_plotly, plot_components_plotly, add_changepoints_to_plot\nfrom sklearn.metrics import mean_squared_error\nfrom statsmodels.tsa.seasonal import seasonal_decompose\nimport statsmodels.api as sm\nimport calendar\n","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:25:29.195166Z","iopub.execute_input":"2022-07-11T18:25:29.195772Z","iopub.status.idle":"2022-07-11T18:25:31.912802Z","shell.execute_reply.started":"2022-07-11T18:25:29.195735Z","shell.execute_reply":"2022-07-11T18:25:31.911787Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.io as pio\npio.templates[\"draft\"] = go.layout.Template(\n    layout_annotations=[\n        dict(\n            textangle=-30,\n            opacity=0.1,\n            font=dict(color=\"black\", size=100),\n            xref=\"paper\",\n            yref=\"paper\",\n            x=0.5,\n            y=0.5,\n            showarrow=False,\n        )\n    ]\n)\npio.templates.default = \"draft\"","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:25:31.914944Z","iopub.execute_input":"2022-07-11T18:25:31.915394Z","iopub.status.idle":"2022-07-11T18:25:31.969295Z","shell.execute_reply.started":"2022-07-11T18:25:31.915347Z","shell.execute_reply":"2022-07-11T18:25:31.968266Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.set_option('display.max_columns', 2000)","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:25:31.970771Z","iopub.execute_input":"2022-07-11T18:25:31.971101Z","iopub.status.idle":"2022-07-11T18:25:31.975413Z","shell.execute_reply.started":"2022-07-11T18:25:31.971068Z","shell.execute_reply":"2022-07-11T18:25:31.974336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","execution":{"iopub.status.busy":"2022-07-11T18:38:01.351099Z","iopub.execute_input":"2022-07-11T18:38:01.351706Z","iopub.status.idle":"2022-07-11T18:38:01.360326Z","shell.execute_reply.started":"2022-07-11T18:38:01.351650Z","shell.execute_reply":"2022-07-11T18:38:01.359281Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- `calendar.csv` - Contains information about the dates on which the products are sold.\n- `sales_train_validation.csv` - Contains the historical daily unit sales data per product and store [d_1 - d_1913]\n- `sample_submission.csv` - The correct format for submissions. Reference the Evaluation tab for more info.\n- `sell_prices.csv` - Contains information about the price of the products sold per store and date.\n- `sales_train_evaluation.csv` - Includes sales [d_1 - d_1941] (labels used for the Public leaderboard)","metadata":{}},{"cell_type":"markdown","source":"# 1. Reading & Processing Data","metadata":{}},{"cell_type":"code","source":"calendar_data = pd.read_csv('/kaggle/input/m5-forecasting-accuracy/calendar.csv')\ncalendar_data.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:38:02.018440Z","iopub.execute_input":"2022-07-11T18:38:02.019140Z","iopub.status.idle":"2022-07-11T18:38:02.047055Z","shell.execute_reply.started":"2022-07-11T18:38:02.019093Z","shell.execute_reply":"2022-07-11T18:38:02.046119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sales_train_validation = pd.read_csv('/kaggle/input/m5-forecasting-accuracy/sales_train_validation.csv')\nsales_train_validation.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:38:02.231566Z","iopub.execute_input":"2022-07-11T18:38:02.231950Z","iopub.status.idle":"2022-07-11T18:38:07.706524Z","shell.execute_reply.started":"2022-07-11T18:38:02.231913Z","shell.execute_reply":"2022-07-11T18:38:07.705431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sales_train_evaluation = pd.read_csv('/kaggle/input/m5-forecasting-accuracy/sales_train_evaluation.csv')\nsales_train_evaluation.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:38:07.708476Z","iopub.execute_input":"2022-07-11T18:38:07.708763Z","iopub.status.idle":"2022-07-11T18:38:12.796730Z","shell.execute_reply.started":"2022-07-11T18:38:07.708735Z","shell.execute_reply":"2022-07-11T18:38:12.795792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Total NaNs Values in Validation Dataset: {sales_train_validation.isna().sum().sum()}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:38:12.798349Z","iopub.execute_input":"2022-07-11T18:38:12.798634Z","iopub.status.idle":"2022-07-11T18:38:12.934631Z","shell.execute_reply.started":"2022-07-11T18:38:12.798606Z","shell.execute_reply":"2022-07-11T18:38:12.933405Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"initial_memory = sales_train_validation.memory_usage().sum()/(1024**2)\nd_cols = [col for col in sales_train_validation.columns if \"d_\" in col]\n### sales_train_validation[d_cols].max().max()\nsales_train_validation[d_cols] = sales_train_validation[d_cols].astype(np.int32)   # int16 may create issues\nfinal_memory = sales_train_validation.memory_usage().sum()/(1024**2)\nmemory_reduction = (initial_memory - final_memory)/ initial_memory\nprint(f\"Memory Reduction: {round(initial_memory - final_memory, 2)}MB ({round(100*memory_reduction, 2)}%)\")","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:38:12.936050Z","iopub.execute_input":"2022-07-11T18:38:12.936344Z","iopub.status.idle":"2022-07-11T18:41:39.115484Z","shell.execute_reply.started":"2022-07-11T18:38:12.936314Z","shell.execute_reply":"2022-07-11T18:41:39.114124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. Exploratory Data Analysis","metadata":{}},{"cell_type":"markdown","source":"<i>`sales_train_validation` is the train dataset -- Days Range :[D1 - D1913].</i> <br>\n<i>`sales_train_evalutaion` is the evaluation dataset -- Days Range :[D1914 - D1941].</i>","metadata":{}},{"cell_type":"code","source":"train_data = sales_train_validation.copy()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:41:39.118361Z","iopub.execute_input":"2022-07-11T18:41:39.118717Z","iopub.status.idle":"2022-07-11T18:41:39.229535Z","shell.execute_reply.started":"2022-07-11T18:41:39.118682Z","shell.execute_reply":"2022-07-11T18:41:39.228685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.1 General Information","metadata":{}},{"cell_type":"code","source":"print(\"In the dataset, there are:\")\nprint(f\"- {train_data.state_id.nunique()} states: {train_data.state_id.unique().tolist()}.\")\nprint(f\"- {train_data.store_id.nunique()} stores: {train_data.store_id.unique().tolist()}.\")\nprint(f\"- {train_data.cat_id.nunique()} categories: {train_data.cat_id.unique()}.\")\nprint(f\"- {train_data.dept_id.nunique()} departments: {train_data.dept_id.unique().tolist()}.\")\nprint(f\"- {train_data.item_id.nunique()} items.\")","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:53:26.436801Z","iopub.execute_input":"2022-07-11T18:53:26.437191Z","iopub.status.idle":"2022-07-11T18:53:26.463237Z","shell.execute_reply.started":"2022-07-11T18:53:26.437153Z","shell.execute_reply":"2022-07-11T18:53:26.462080Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Total Numner of Time Series: {train_data.id.nunique()}.\")","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:54:06.329713Z","iopub.execute_input":"2022-07-11T18:54:06.330080Z","iopub.status.idle":"2022-07-11T18:54:06.345062Z","shell.execute_reply.started":"2022-07-11T18:54:06.330050Z","shell.execute_reply":"2022-07-11T18:54:06.344208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train_data[:10]","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:54:50.285835Z","iopub.execute_input":"2022-07-11T18:54:50.286436Z","iopub.status.idle":"2022-07-11T18:54:50.291280Z","shell.execute_reply.started":"2022-07-11T18:54:50.286399Z","shell.execute_reply":"2022-07-11T18:54:50.290142Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.2 Example of Time Series","metadata":{}},{"cell_type":"code","source":"non_d_cols = list(set(sales_train_validation.columns) - set(d_cols))","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:55:34.938587Z","iopub.execute_input":"2022-07-11T18:55:34.938962Z","iopub.status.idle":"2022-07-11T18:55:34.943804Z","shell.execute_reply.started":"2022-07-11T18:55:34.938925Z","shell.execute_reply":"2022-07-11T18:55:34.942954Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_data.loc[train_data.d_100 == train_data.d_100.max()]","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:55:45.356237Z","iopub.execute_input":"2022-07-11T18:55:45.356629Z","iopub.status.idle":"2022-07-11T18:55:46.122965Z","shell.execute_reply.started":"2022-07-11T18:55:45.356592Z","shell.execute_reply":"2022-07-11T18:55:46.122104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's take the example of item <b>FOODS_3_586</b> sales in California store <b>TX_3</b>","metadata":{}},{"cell_type":"code","source":"def merge_with_calendar(data, calendar_data):\n    # data should have a date column \"d\"\n    assert 'd' in data.columns, 'DataFrame should have a column \"d\" !'\n    # Merge With Calendar\n    cal = calendar_data[['d', 'date', 'wm_yr_wk', 'weekday', 'wday', 'month', 'year',\n                    'event_name_1', 'event_type_1', 'event_name_2', 'event_type_2']]\n    d = pd.merge(cal, data, on=\"d\")\n    # Fill Missing Event Values with None\n    for col in ['event_name_1', 'event_type_1', 'event_name_2', 'event_type_2']:\n        d[col].fillna('None', inplace=True)\n    return d","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:56:25.062839Z","iopub.execute_input":"2022-07-11T18:56:25.063500Z","iopub.status.idle":"2022-07-11T18:56:25.071423Z","shell.execute_reply.started":"2022-07-11T18:56:25.063442Z","shell.execute_reply":"2022-07-11T18:56:25.070348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_ts_example(data, d_cols, calendar_data, item_id, store_id, idx=None):\n    try:\n        if idx is None:\n            ts = data.loc[(data['item_id'] == item_id) & (data['store_id'] == store_id)]\n            ts = ts[d_cols].T.reset_index()\n            ts.columns = ['d', 'sales']\n        else:\n            ts = data.loc[idx][d_cols].reset_index()\n            ts.columns = ['d', 'sales']\n        # Make sure that sales column's type is int\n        ts[\"sales\"] = ts[\"sales\"].astype(\"int\")\n        return merge_with_calendar(ts, calendar_data)\n    except Exception as e:\n        print(f'Can not extract time series: {e}')","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:56:27.616715Z","iopub.execute_input":"2022-07-11T18:56:27.617152Z","iopub.status.idle":"2022-07-11T18:56:27.625499Z","shell.execute_reply.started":"2022-07-11T18:56:27.617113Z","shell.execute_reply":"2022-07-11T18:56:27.624369Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ts = get_ts_example(train_data, d_cols, calendar_data, item_id='FOODS_3_586', store_id='TX_3')","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:56:36.358929Z","iopub.execute_input":"2022-07-11T18:56:36.359513Z","iopub.status.idle":"2022-07-11T18:56:36.393192Z","shell.execute_reply.started":"2022-07-11T18:56:36.359457Z","shell.execute_reply":"2022-07-11T18:56:36.392151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ts","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:56:36.747787Z","iopub.execute_input":"2022-07-11T18:56:36.748211Z","iopub.status.idle":"2022-07-11T18:56:36.773032Z","shell.execute_reply.started":"2022-07-11T18:56:36.748178Z","shell.execute_reply":"2022-07-11T18:56:36.772054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Length of the time series: {len(ts)} days.\")\nprint(f\"Years Available in the time series: {ts.year.unique()}\")","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:57:27.938559Z","iopub.execute_input":"2022-07-11T18:57:27.938946Z","iopub.status.idle":"2022-07-11T18:57:27.945227Z","shell.execute_reply.started":"2022-07-11T18:57:27.938909Z","shell.execute_reply":"2022-07-11T18:57:27.944394Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(ts['event_type_2'].value_counts())","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:57:32.567613Z","iopub.execute_input":"2022-07-11T18:57:32.567996Z","iopub.status.idle":"2022-07-11T18:57:32.575816Z","shell.execute_reply.started":"2022-07-11T18:57:32.567960Z","shell.execute_reply":"2022-07-11T18:57:32.574899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<i>From 1913 days there are only 4 days with type 2 event! As a result we'll ignore that event type.</i>","metadata":{}},{"cell_type":"markdown","source":"## 2.3 Plots","metadata":{"execution":{"iopub.status.busy":"2022-07-11T19:00:00.154422Z","iopub.execute_input":"2022-07-11T19:00:00.155024Z","iopub.status.idle":"2022-07-11T19:00:00.368094Z","shell.execute_reply.started":"2022-07-11T19:00:00.154978Z","shell.execute_reply":"2022-07-11T19:00:00.366699Z"}}},{"cell_type":"code","source":"events_1_data = ts['event_type_1'].value_counts().iloc[1:]","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:59:49.989470Z","iopub.execute_input":"2022-07-11T18:59:49.989800Z","iopub.status.idle":"2022-07-11T18:59:49.998426Z","shell.execute_reply.started":"2022-07-11T18:59:49.989770Z","shell.execute_reply":"2022-07-11T18:59:49.997370Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(11, 5))\nax = fig.add_subplot()\nax.pie(x=events_1_data.values,\n       labels=events_1_data.index,\n       shadow=True,\n       radius=1,\n       autopct='%1.1f%%')\nax.set_title('Distribution of Type 1 Events')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T18:58:43.253253Z","iopub.execute_input":"2022-07-11T18:58:43.253683Z","iopub.status.idle":"2022-07-11T18:58:43.337740Z","shell.execute_reply.started":"2022-07-11T18:58:43.253644Z","shell.execute_reply":"2022-07-11T18:58:43.336846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Alternative WIth Plotly\n# Target Variable\nfig = go.Figure()\n\nfig.add_trace(go.Pie(labels=events_1_data.index,\n                     values=events_1_data.values,\n                     textfont_size=20, marker=dict(line=dict(color='#000000', width=2)))\n             )\n\nfig.update_layout(title=\"Events Type 1 Distribution\")\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T19:01:43.580310Z","iopub.execute_input":"2022-07-11T19:01:43.580731Z","iopub.status.idle":"2022-07-11T19:01:43.788750Z","shell.execute_reply.started":"2022-07-11T19:01:43.580693Z","shell.execute_reply":"2022-07-11T19:01:43.787639Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(20, 6))\nax = fig.add_subplot()\nax.plot(ts.date, ts.sales)\nax.set_xticks(ts.date.values[::90])\nax.set_xlabel('Date')\nax.set_ylabel('Sales')\nax.grid()\nax.set_title(f'Daily Sales of FOODS_3_586 in TX_3 Store')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T19:05:25.465937Z","iopub.execute_input":"2022-07-11T19:05:25.466325Z","iopub.status.idle":"2022-07-11T19:05:25.931116Z","shell.execute_reply.started":"2022-07-11T19:05:25.466295Z","shell.execute_reply":"2022-07-11T19:05:25.929920Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Alternative ith Plotly\nfig = go.Figure()\nfig.add_trace(go.Scatter(x=ts.date, y=ts.sales))\nfig.update_layout(title=\"Daily Sales of FOODS_3_586 in TX_3 Store\", xaxis_title=\"Date\", yaxis_title=\"Sales\")\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T19:04:09.536827Z","iopub.execute_input":"2022-07-11T19:04:09.537230Z","iopub.status.idle":"2022-07-11T19:04:09.575230Z","shell.execute_reply.started":"2022-07-11T19:04:09.537195Z","shell.execute_reply":"2022-07-11T19:04:09.574013Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Area Plot\nfig = px.area(ts, \n              x='date',\n              y='sales',\n              title='Time Series Yearly Area Plot',\n              facet_row='year',\n              facet_row_spacing=0.05)\nfig.update_layout(width=900,\n                 height=900)\n             \nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T19:05:35.728584Z","iopub.execute_input":"2022-07-11T19:05:35.728970Z","iopub.status.idle":"2022-07-11T19:05:36.051778Z","shell.execute_reply.started":"2022-07-11T19:05:35.728932Z","shell.execute_reply":"2022-07-11T19:05:36.050744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for year in ts.year.unique():\n    t = ts[ts.year == year]\n    fig = plt.figure(figsize=(20, 6))\n    ax = fig.add_subplot()\n    t_event_1 = t.loc[t.event_type_1 != 'None']\n    ax.plot(t.date, t.sales)\n    ax.scatter(t_event_1.date,\n               t_event_1.sales,\n               color='red',\n               label='Type 1 Event')\n    ax.set_xticks(t.date.values[::30])\n    ax.set_xlabel('Date')\n    ax.set_ylabel('Sales')\n    ax.grid()\n    ax.set_title(f'Sales of FOODS_3_586 in TX_3 store for {year}')\n    ax.legend()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T19:06:01.238117Z","iopub.execute_input":"2022-07-11T19:06:01.238471Z","iopub.status.idle":"2022-07-11T19:06:02.700541Z","shell.execute_reply.started":"2022-07-11T19:06:01.238441Z","shell.execute_reply":"2022-07-11T19:06:02.699676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(20, 15))\nax1 = fig.add_subplot(311)\nsns.boxplot(data=ts, x='month', y='sales', ax=ax1)\nplt.grid(axis=\"y\")\nax2 = fig.add_subplot(312)\nsns.boxplot(data=ts, x='weekday', y='sales', ax=ax2)\nplt.grid(axis=\"y\")\nax3 = fig.add_subplot(313)\nsns.boxplot(data=ts, x='event_type_1', y='sales', ax=ax3)\nplt.grid(axis=\"y\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T19:07:18.492172Z","iopub.execute_input":"2022-07-11T19:07:18.492686Z","iopub.status.idle":"2022-07-11T19:07:19.426305Z","shell.execute_reply.started":"2022-07-11T19:07:18.492651Z","shell.execute_reply":"2022-07-11T19:07:19.425460Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- A simple observation: from the first plot we can see that the biggest part of sales of this item (<i>FOODS_3_586</i>) in this store (<i>TX_3</i>) is during Month 8 (August).<br>\n- This very detailed level (<i>Item level</i>) won't generate many insights, aggregated levels will do such as <i>State</i>, <i>Store</i>, <i>Category</i> and <i>Department</i> levels.","metadata":{}},{"cell_type":"markdown","source":"Plotly offers some interesting plots and visualizations ! Let's try some of them at our time series.","metadata":{}},{"cell_type":"code","source":"fig = px.histogram(ts,\n                   x='sales',\n                   marginal='box',\n                   title='Sales Distribution for FOODS_3_586 at TX_3 store')\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T19:07:48.890910Z","iopub.execute_input":"2022-07-11T19:07:48.891287Z","iopub.status.idle":"2022-07-11T19:07:49.021219Z","shell.execute_reply.started":"2022-07-11T19:07:48.891250Z","shell.execute_reply":"2022-07-11T19:07:49.020148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.line(data_frame=ts, x='date',  y='sales', \n              color='year', \n              title='Sales of FOODS_3_586 at TX_3 Store',)\nfig.update_layout(legend=dict(x=0.25, y=-0.2, title_font_family=\"Times New Roman\", orientation=\"h\",\n                              bgcolor=\"snow\", bordercolor=\"Black\",borderwidth=1),font=dict(size=11)\n                 )\n# plot and legend are Interactive !\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T19:09:44.128778Z","iopub.execute_input":"2022-07-11T19:09:44.129176Z","iopub.status.idle":"2022-07-11T19:09:44.212130Z","shell.execute_reply.started":"2022-07-11T19:09:44.129140Z","shell.execute_reply":"2022-07-11T19:09:44.210896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"An interesting pattern to observe is that, for this item example, sales decrease to the lowest level at the end of each year !","metadata":{}},{"cell_type":"code","source":"for feature in ['wday', 'month', 'year', 'event_name_1']: #wday = weekday\n    data_feature = ts.groupby(feature).mean()['sales'].reset_index().sort_values(by='sales')\n    fig = px.bar(data_frame=data_feature,\n          x=feature,\n          y='sales',\n          title=f'Average Sales by {feature}')\n    fig.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-11T19:10:01.053005Z","iopub.execute_input":"2022-07-11T19:10:01.053345Z","iopub.status.idle":"2022-07-11T19:10:01.240916Z","shell.execute_reply.started":"2022-07-11T19:10:01.053316Z","shell.execute_reply":"2022-07-11T19:10:01.239953Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can easily notice that the item's Average sale is at its maximum at Father's Day.","metadata":{}},{"cell_type":"markdown","source":"## 2.4 Aggregated Level Analysis","metadata":{}},{"cell_type":"markdown","source":"- Now we'll go through more aggregated analysis.","metadata":{}},{"cell_type":"markdown","source":" - The Data hierarchy is:<br>\n <b>State</b> ==>  <b>Store</b> ==>  <b>Category</b> ==>  <b>Department</b> ==>  <b>Item</b>","metadata":{}},{"cell_type":"code","source":"d = train_data.groupby(['state_id','store_id']).count().reset_index()\nd['state_id'].value_counts().plot(kind=\"bar\", grid=True, title=\"(A) Number of Stores by State\", yticks=[0,1,2,3,4])\nplt.show()\nd = train_data.groupby(['state_id','id']).count().reset_index()\nd['state_id'].value_counts().plot(kind=\"bar\", grid=True, title=\"(B) Number of Items by State\")\nplt.show()\nd = train_data.groupby(['store_id','item_id']).count().reset_index()\nd['store_id'].value_counts().plot(kind=\"bar\", grid=True, title=\" (C) Number of Items by Store\")\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- We will use <b>pandas.DataFrame.stack</b> function to transform our data.\n- The result data should have as columns <b>['d', 'state_id','store_id','cat_id','dept_id','item_id', 'id', 'sales']</b> where \"d\" column has values in <br>[d_1,..., d_1913].\n- Initially we have data with 10 stores and 3049 items per store so as a result we have 30490 time series !\n- Each time series contains sales data of 1913 days. Final output data will have then 8 columns and 30490*1913 = 58327370 rows !","metadata":{}},{"cell_type":"code","source":"train_data.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Transform Data Structure\ndata = train_data.set_index(non_d_cols)\n# the following will make one column for sales and one columns for \"d\" values (d_1 ... d_1913)\ndata = data.stack()\ndata = data.to_frame() \ndata.columns = [\"sales\"]\ndata.reset_index(inplace=True)\ndata.columns = non_d_cols + [\"d\", \"sales\"]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<i>The followoing are function that will be used for different analysis.</i>","metadata":{}},{"cell_type":"code","source":"def plot_daily_data(df,level):\n    levels_dict = {'cat':'Category', 'dept':'Department', 'store':'Store', 'state':'State'}\n    fig = px.line(data_frame=df, \n                  x='date', \n                  y='sales', \n                  color=f'{level}_id', \n                  title=f'Sales by {levels_dict[level]}')\n    fig.update_layout(legend=dict(x=1,\n                                  y=1,\n                                  title_font_family=\"Times New Roman\",\n                                  bgcolor=\"snow\",\n                                  bordercolor=\"Black\",\n                                  borderwidth=1),\n                      font=dict(size=11)\n                     )\n    fig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_yearly_data(df,level):\n    levels_dict = {'cat':'Category', 'dept':'Department', 'store':'Store', 'state':'State'}\n    n_years = df.year.nunique()\n    years = df.year.unique()\n    level_elements = df[f'{level}_id'].unique().tolist()  \n    all_colors= ['red','green','blue','purple','cyan','orange','pink','yellow','black','magenta']\n    colors = all_colors[:len(level_elements)]\n    fig = plt.figure(figsize=(15,15))\n    gs = GridSpec(n_years, len(level_elements))\n    c_idx = 0\n    for l_e, color in zip(level_elements,colors):\n        r_idx = 0\n        for year in years:\n            ax = fig.add_subplot(gs[r_idx,c_idx], xticks=[], yticks=[])\n            df1 = df.loc[(df.year == year) & (df[f'{level}_id'] == l_e)]\n            ax.plot(df1.date, df1['sales'], color=color, linewidth=0.9)\n            ax.set_title(f'{l_e}: {year}')\n            r_idx += 1\n        c_idx+=1\n    fig.suptitle(f'Yearly Sales by {levels_dict[level]}')\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_average_sales(df,level):\n    for feature in ['weekday', 'month', 'year', 'event_name_1']:\n        data_feature = df.groupby([f'{level}_id', feature]).mean()['sales'].reset_index()\n        fig = px.bar(data_frame=data_feature,\n              x=feature,\n              y='sales',\n              color=f'{level}_id',\n              title=f'Average Sales by {feature}')\n        fig .update_layout(legend=dict(x=1,\n                                       y=1,\n                                       title_font_family=\"Times New Roman\",\n                                       bgcolor=\"mintcream\",\n                                       bordercolor=\"black\",\n                                       borderwidth=1))\n        fig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<b>For plotly plots, you can double click on legend to visualize data parts separately. You can also zoom in, zoom out and autoscale plots.</b>","metadata":{}},{"cell_type":"markdown","source":"### 2.4.1 State Level Analysis","metadata":{}},{"cell_type":"code","source":"# Preparing state-level data\nstate_data = data.groupby([\"state_id\",\"d\"]).sum()[\"sales\"]\nstate_data = state_data.reset_index()\nstate_data = merge_with_calendar(state_data, calendar_data)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"state_data.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_daily_data(state_data,'state')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.histogram(state_data,\n                   x='sales',\n                   color='state_id',\n                   marginal='box',\n                   title='Sales Distribution By State')\nfig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- California has the highest number of sales.\n- We observe again the sales decrease to their lowest level (previously observed with a single item time series) at the end of every year, let's try to get more details about this pattern !","metadata":{}},{"cell_type":"code","source":"dlow = state_data.loc[(state_data['sales']<20) & (state_data['month']==12)]\ndlow.style.applymap(lambda x:\"background-color:yellow\", subset=['event_name_1'])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<b> It was Christmas effect !</b>","metadata":{}},{"cell_type":"code","source":"plot_yearly_data(state_data,'state')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_average_sales(state_data,'state')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Sales average is increasing over years.\n- Weekends correspond to the highest sales average.\n- Sales average is almost the same over different months.","metadata":{}},{"cell_type":"markdown","source":"### 2.4.2 Store Level Analysis","metadata":{}},{"cell_type":"code","source":"# Preparing store-level data\nstore_data = data.groupby([\"store_id\",\"d\"]).sum()[\"sales\"]\nstore_data = store_data.reset_index()\nstore_data = merge_with_calendar(store_data, calendar_data)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"store_data.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_daily_data(store_data,'store')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.histogram(store_data,\n                   x='sales',\n                   color='store_id',\n                   marginal='box',\n                   title='Sales Distribution By Store')\nfig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- CA_3 is the store having the highest number of sales.","metadata":{}},{"cell_type":"code","source":"plot_yearly_data(store_data,'store')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_average_sales(store_data,'store')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Let's see Sales' correlations between different Stores.","metadata":{}},{"cell_type":"code","source":"corr_data = pd.pivot_table(data=store_data,\n                           index='date',\n                           values='sales',\n                           columns='store_id')\ncorr_data.sort_values(by=\"date\", inplace=True)\nplt.figure(figsize=(12,5))\nheatmap = sns.heatmap(corr_data.corr(), annot=True, fmt='.2f')\nheatmap.set_yticklabels(heatmap.yaxis.get_ticklabels(), rotation=0)\nheatmap.set_title('Correlation Between Stores Sales')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- The highest correlation is between CA_1 and CA_3 stores, in the same state.\n- The lowest correlation is between WI_1 and WI_3 stores, in the same state!","metadata":{}},{"cell_type":"markdown","source":"### 2.4.3 Catgeory Level Analysis","metadata":{}},{"cell_type":"code","source":"train_data.groupby('cat_id').count()['id'].reset_index().plot(x='cat_id', \n                                                              kind='bar', \n                                                              figsize=(15,5),\n                                                              grid=True,\n                                                              title='Number of Items by Category')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Preparing category-level data\ncat_data = data.groupby([\"cat_id\",\"d\"]).sum()[\"sales\"]\ncat_data = cat_data.reset_index()\ncat_data = merge_with_calendar(cat_data, calendar_data)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat_data.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_daily_data(cat_data,'cat')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- FOODS is the category having the highest number of sales, HOBBIES having the lowest one.","metadata":{}},{"cell_type":"code","source":"fig = px.histogram(cat_data,\n                   x='sales',\n                   color='cat_id',\n                   marginal='box',\n                   title='Sales Distribution By Category')\nfig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_yearly_data(cat_data,'cat')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_average_sales(cat_data,'cat')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's plot a Sales Calendar Heatmap of the first and last years (by Categroy).","metadata":{}},{"cell_type":"code","source":"years = cat_data.year.unique()\nfor year in [years[0],years[-1]]:\n    for cat in cat_data.cat_id.unique():\n        dyear = cat_data.loc[(cat_data.year == year) & (cat_data.cat_id == cat)] # & (cat_data.cat_id == cat)\n        fig = go.Figure(data=go.Heatmap(\n                z=dyear.sales, #z,\n                x=dyear.date, #dates,\n                y=dyear.cat_id, #programmers,\n                colorscale=px.colors.sequential.Plasma_r))\n\n        fig.update_layout(\n            title=f'{cat} Sales {year}',\n            xaxis_nticks=36)\n\n        fig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.4.4 Department Level Analysis","metadata":{}},{"cell_type":"code","source":"train_data.groupby('dept_id').count()['id'].reset_index().plot(x='dept_id', \n                                                              kind='bar', \n                                                              figsize=(15,5),\n                                                              grid=True,\n                                                              title='Number of Items by Department')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Preparing department-level data\ndept_data = data.groupby([\"dept_id\",\"d\"]).sum()[\"sales\"]\ndept_data = dept_data.reset_index()\ndept_data = merge_with_calendar(dept_data, calendar_data)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dept_data.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_daily_data(dept_data,'dept')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- FOODS_3 and HOBBIES_2 have respectively the highest and lowest number of sales among all departments.","metadata":{}},{"cell_type":"code","source":"fig = px.histogram(dept_data,\n                   x='sales',\n                   color='dept_id',\n                   marginal='box',\n                   title='Sales Distribution By Department')\nfig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_yearly_data(dept_data,'dept')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_average_sales(dept_data,'dept')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"corr_data = pd.pivot_table(data=dept_data,\n                           index='date',\n                           values='sales',\n                           columns='dept_id')\ncorr_data.sort_values(by=\"date\", inplace=True)\nplt.figure(figsize=(12,5))\nheatmap = sns.heatmap(corr_data.corr(), annot=True, fmt='.2f')\nheatmap.set_yticklabels(heatmap.yaxis.get_ticklabels(), rotation=0)\nheatmap.set_title('Correlation Between Department Sales')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"After the analyzing and visualizing part comes the forecast part  !","metadata":{}},{"cell_type":"markdown","source":"# 3. Forecast","metadata":{}},{"cell_type":"markdown","source":"- In this part, we'll choose 6 Time Series from the 30490 ones we have and use them to test different forecasting approaches and evaluate them.","metadata":{}},{"cell_type":"code","source":"forecast_horizon = 28 # from d_1914 to d_1941\nd_fcst_columns = sales_train_evaluation.columns[-forecast_horizon:].tolist()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_ground_truth(idx, df, d_fcst_columns):\n    return df.loc[idx, d_fcst_columns].values  ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_results(fcst, y_eval, rmse, algo, item):\n    fig = plt.figure(figsize=(11, 5))\n    ax = fig.add_subplot()\n    ax.plot(fcst, color='red', label='Forecast')\n    ax.plot(y_eval, color='blue', label='Ground Truth')\n    ax.set_title(f' {algo} for {item}, RMSE: {rmse}')\n    ax.grid()\n    ax.legend()\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.1 Time Series Selection","metadata":{}},{"cell_type":"markdown","source":"We will choose:\n- 3 Time Series with enough data and few zeros.\n- 3 Time Series with many zeros.","metadata":{}},{"cell_type":"code","source":"choice_data = train_data.copy()\nchoice_data['d_val'] = choice_data[d_cols].mean(axis=1)\nchoice_data.drop(columns=d_cols,inplace=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f'Time Series Averages: Min {choice_data.d_val.min()} Max:{choice_data.d_val.max()} Median: {choice_data.d_val.median()}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"idx1 = choice_data.loc[choice_data['d_val'] >= 50].sample(n=3, random_state=1).index\nidx2 = choice_data.loc[(choice_data['d_val'] <= 5) & (choice_data['d_val'] > 1)].sample(n=3, random_state=1).index\nidx = idx1.tolist() + idx2.tolist()\nts_test = train_data.iloc[idx]#.reset_index(drop=True)\ntest_items = ts_test.id.unique().tolist()\ntest_items = [x[:-11] for x in test_items]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_items","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Dataframe for RMSE\nrmse_summary = pd.DataFrame({\"items\":test_items}, index=idx)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.2 Statistical Method: ARIMA","metadata":{}},{"cell_type":"markdown","source":"<b>ARIMA</b> stands for <b>A</b>uto<b>R</b>egressive <b>I</b>ntegrated <b>M</b>oving <b>A</b>verage.","metadata":{}},{"cell_type":"markdown","source":"Types:\n- <b>ARIMA</b>: Non-Seasonal. \n- <b>SARIMA</b>: Seasonal ARIMA.\n- <b>SARIMAX</b>: Seasonal ARIMA with eXogenous variables.","metadata":{}},{"cell_type":"code","source":"def choose_sarimax_order_and_forecast(ps,ds,qs, y_train, y_eval):\n    best_model, best_rmse, best_order, best_fcst = None, None, None, None\n    for p in ps:\n        for d in ds:\n            for q in qs:\n                order = (p,d,q)\n                model = sm.tsa.SARIMAX(y_train, \n                               order=order, \n                               trend='c',\n                               enforce_invertibility=False,\n                               enforce_stationarity=False).fit(disp=False, warn_convergence=False)\n                fcst = model.predict(start=len(y_train), end=len(y_train) - 1 + len(y_eval))\n                try:\n                    fcst = [round(x) for x in fcst]\n                    rmse = round(np.sqrt(mean_squared_error(fcst, y_eval)), 3)\n                    if (best_rmse is None) or (rmse < best_rmse):\n                        best_model, best_rmse, best_order, best_forecast= model, rmse, order, fcst\n                except Exception as e:\n                    print(f'For order={order}, model results are invalid: {e}')\n    print(f\"Best Order: {best_order}\")\n    return best_rmse, best_forecast","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df0 = get_ts_example(ts_test, d_cols, calendar_data, item_id=None, store_id=None, idx=idx[0])\ndf0['date'] = df0['date'].apply(lambda x : pd.to_datetime(x))\ndf0 = df0[['date','sales']]\ny_train = df0[\"sales\"].values\ndf0","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"result = seasonal_decompose(y_train, model='additive', period=365)\nfig = result.plot()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(2,1,figsize=(20,7))\nsm.tsa.graphics.plot_acf(y_train, lags=30, ax=ax[0])\nax[0].set_title('Autocorreation Function: lags=30')\nsm.tsa.graphics.plot_pacf(y_train, lags=30, ax=ax[1])\nax[1].set_title('Partial Autocorreation Function: lags=30')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"order = (p,d,q)\n- p: AutoRegression (AR) order.\n- d: Trend Differncing order.\n- q: Moving Average (MA) order. <br>","metadata":{}},{"cell_type":"code","source":"# Plotting Autocorrelation with pandas\nfig, ax = plt.subplots(1,1,figsize=(20,7))\npd.plotting.autocorrelation_plot(y_train, ax=ax)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sarimax_model = sm.tsa.SARIMAX(y_train, \n                               order=(7,1,7), \n                               trend='c',\n                               enforce_invertibility=False,\n                               enforce_stationarity=False).fit(disp=False, warn_convergence=False)\nsarimax_model.summary()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_eval = get_ground_truth(idx[0], sales_train_evaluation, d_fcst_columns)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fcst = sarimax_model.predict(start=len(y_train), end=len(y_train) - 1 + len(y_eval))\nrmse = round(np.sqrt(mean_squared_error(fcst, y_eval)), 3)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_results(fcst, y_eval, rmse, \"SARIMAX\", test_items[0])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# SARIMAX Parameters Grid\nps = range(1,8)\nds = range(0,2)\nqs = range(0,8)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rmse_sarimax = []\nfor i,ix in enumerate(idx):\n    print(f\"Processing {test_items[i]}...\")\n    # Get Time Series (Train)\n    df0 = get_ts_example(ts_test, d_cols, calendar_data, item_id=None, store_id=None, idx=ix)\n    df0['date'] = df0['date'].apply(lambda x : pd.to_datetime(x))\n    df0 = df0[['date','sales']]\n    y_train = df0[\"sales\"].values\n    y_eval = get_ground_truth(ix, sales_train_evaluation, d_fcst_columns)\n    # Train SARIMAX model\n    rmse, fcst = choose_sarimax_order_and_forecast(ps,ds,qs, y_train, y_eval)\n    # Plot\n    plot_results(fcst, y_eval, rmse, \"SARIMAX\", test_items[i])\n    rmse_sarimax.append(rmse)","metadata":{"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Dataframe for RMSE\nrmse_summary = pd.DataFrame({\"items\":test_items}, index=idx)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rmse_summary[\"RMSE_SARIMAX\"] = rmse_sarimax\nrmse_summary","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.3 FB Prophet","metadata":{}},{"cell_type":"markdown","source":"FB Prophet requires a column 'ds' (for dates) and a columns 'y' (target variable).","metadata":{}},{"cell_type":"code","source":"def generate_holidays(calendar_data, dates):\n    holidays = calendar_data.loc[calendar_data['d'].isin(dates), ['date','event_name_1']].dropna()\n    holidays['ds'] = holidays['date'].apply(lambda x : pd.to_datetime(x))\n    holidays['upper_window'] = 0\n    holidays['lower_window'] = 0\n    holidays.rename(columns={\"event_name_1\":\"holiday\"}, inplace=True)\n    holidays.drop(columns='date', inplace=True)\n    holidays.reset_index(drop=True, inplace=True)\n    return holidays","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# holidays for FB prophet model\nholidays = generate_holidays(calendar_data, d_cols+d_fcst_columns)\nholidays","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = get_ts_example(ts_test, d_cols, calendar_data, item_id=None, store_id=None, idx=idx[0])\ndf['ds'] = df['date'].apply(lambda x : pd.to_datetime(x))\ndf = df[['ds','sales']]\ndf = df.rename(columns={'sales':'y'})\ndf","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['cap'] = df['y'].max()\ndf['floor'] = df['y'].min()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Creating Model\nmodel = Prophet(daily_seasonality=True, \n                holidays=holidays,\n                holidays_prior_scale=0.2,\n                growth='logistic', # possible value: 'logistic', 'linear' or 'flat'\n                changepoint_range=0.8, # default value\n                changepoint_prior_scale=0.2) # default value\nmodel.add_seasonality(name='weekly', period=7, fourier_order=3)\nmodel.add_seasonality(name='yearly', period=364, fourier_order=10)\nmodel.fit(df)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Forecasting\nfuture = model.make_future_dataframe(periods=forecast_horizon)\nfuture['cap'] = df['y'].max()\nfuture['floor'] = df['y'].min()\nforecast = model.predict(future)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig1 = model.plot(forecast)\na = add_changepoints_to_plot(fig1.gca(), model, forecast)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Interactive plots\nplot_plotly(model, forecast)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Model Components\nplot_components_plotly(model, forecast)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fcst = forecast['yhat'].values[-forecast_horizon:]\nfcst = [int(x) for x in fcst]\nground_truth = get_ground_truth(idx[0], sales_train_evaluation, d_fcst_columns)\nrmse = round(np.sqrt(mean_squared_error(ground_truth, fcst)), 3)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(11, 5))\nax = fig.add_subplot()\nax.plot(fcst, color='red', label='Forecast')\nax.plot(ground_truth, color='blue', label='Ground Truth')\nax.set_title(f'FB Prophet for {test_items[i]}, RMSE: {rmse}')\nax.grid()\nax.legend()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# FB Prophet Parameters Grid\nchangepoint_prior_scale = [0.0001, 0.001, 0.1, 0.5]\nseasonality_prior_scale = [0.01, 0.1, 1, 10]\nholiday_prior_scale = [0.1, 0.2, 0.5, 1]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def tune_fbprophet_and_forecast(y_eval, df, forecast_horizon, chngp, sps, hps):\n    best_rmse, best_fcst = None, None\n    # Choosing the Best Parameters\n    for p1 in chngp:\n        for p2 in sps:\n            for p3 in hps:\n                # Training Model\n                model = Prophet(daily_seasonality=True, \n                                holidays=holidays,\n                                holidays_prior_scale=p3,\n                                seasonality_prior_scale=p2,\n                                growth='logistic', \n                                changepoint_range=0.8, \n                                changepoint_prior_scale=p1)\n                model.add_seasonality(name='weekly', period=7, fourier_order=3)\n                model.add_seasonality(name='yearly', period=364, fourier_order=10)\n                model.fit(df)\n                # Forecasting\n                future = model.make_future_dataframe(periods=forecast_horizon)\n                future['cap'] = df['y'].max()\n                future['floor'] = df['y'].min()\n                forecast = model.predict(future)\n                # Evaluating\n                fcst = forecast['yhat'].values[-forecast_horizon:]\n                fcst = [int(x) for x in fcst]\n                rmse = round(np.sqrt(mean_squared_error(y_eval, fcst)), 3)\n                if (best_rmse is None) or (rmse < best_rmse):\n                    best_rmse, best_fcst = rmse, fcst\n    return best_rmse, best_fcst","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rmse_fbprophet = []\nfor i,ix in enumerate(idx):\n    print(f\"Processing {test_items[i]}...\")\n    # Get Time Series (Train)\n    df = get_ts_example(ts_test, d_cols, calendar_data, item_id=None, store_id=None, idx=ix)\n    df['ds'] = df['date'].apply(lambda x : pd.to_datetime(x))\n    df = df[['ds','sales']]\n    df = df.rename(columns={'sales':'y'})\n    # Add Cap and Floor Columns\n    df['cap'] = df['y'].max()\n    df['floor'] = df['y'].min()\n    # Get Ground Truth\n    y_eval = get_ground_truth(ix, sales_train_evaluation, d_fcst_columns)\n    # Train FBProphet model\n    rmse, fcst = tune_fbprophet_and_forecast(y_eval, df, forecast_horizon, changepoint_prior_scale, seasonality_prior_scale, holiday_prior_scale)\n    # Plot\n    plot_results(fcst, y_eval, rmse, \"FB Prophet\", test_items[i])\n    rmse_fbprophet.append(rmse)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rmse_summary[\"RMSE_FBPROPHET\"] = rmse_fbprophet\nrmse_summary","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"r1 = round(rmse_summary.RMSE_FBPROPHET.mean(), 3)\nr2 = round(rmse_summary.RMSE_SARIMAX.mean(), 3)\nprint(f\"MEAN RMSE SARIMAX: {r2}, FBProphet: {r1}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}