{"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":"# **<span style = 'color:blue'>Neural Basis Expansion Analysis Time Series (NBEATS) Interpretable model for Store sales time series forecasting</span>**\n\n## **<span style='color:green'> Contents:</span>**<a id=\"Table\"></a>\n\n* [Import the libraries](#Import)\n* [Dataset Information](#Dataset)\n* [Visualize the Time Series to study the original plot](#Visualize)\n    > * [Sale Categories of Store Number 1](#Store1)\n    > * [Total Sales Per Year](#Total)\n    > * [Daily Store Sales Data](#Daily)\n* [Feature Engineering & Data Preprocessing](#Preprocessing)\n* [Define Interpretable NBeats model](#Build)\n    > * [Code to structure NBeats architecture](#Architecture)\n    > * [Define NBEATS wrapper to create generic model and interpretable model](#NBEATS-wrapper)\n    > * [Define Generic Block](#Generic)\n    > * [Define Trend Block](#Trend)\n    > * [Define Seasonality Block](#Seasonality)\n    > * [Define Interpretable Block](#Interpretable)\n    > * [Define quantile loss error function](#Loss)\n    > * [Define plot implementation function](#Plot)\n* [Interpretable NBeats for Total Sales](#Total-sale)\n    > * [Model Building & Hypertuning of model](#Building)\n    > * [Model Training](#Training)\n    > * [Output Decomposition](#Decomposition)\n    > * [Model Forecast](#Forecast)\n* [Interpretable NBeats for Daily Store Sales Data](#Store-sale)\n    > * [Model Building & Training](#Building2)\n    > * [Model Forecast of Daily Store Sales](#Forecast2)\n    ","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"# Import data handling & numerical libraries\nimport pandas as pd\nimport numpy as np\nfrom copy import copy\nimport datetime\n\n# Import Data Visualization libraries\nimport seaborn as sb\nimport matplotlib.pyplot as plt\n\n#import libraries for muting unnecessary warnings if needed\nimport warnings\nwarnings.filterwarnings('ignore')","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:05.116945Z","iopub.execute_input":"2022-07-27T08:03:05.117296Z","iopub.status.idle":"2022-07-27T08:03:05.123558Z","shell.execute_reply.started":"2022-07-27T08:03:05.117263Z","shell.execute_reply":"2022-07-27T08:03:05.122251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **<span style = 'color:green'>2. Dataset information</span>**<a id ='Dataset'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\n\nThe aim of this competition is to predict sales for the thousands of product families sold at Favorita stores located in Ecuador. The training data includes dates, store and product information, whether that item was being promoted, as well as the sales numbers. Additional files include supplementary information that may be useful in building various models. \n1. **train.csv**\nThe training data, comprising time series of features store_nbr, family, and onpromotion as well as the target sales.\n    * store_nbr: identifies the store at which the products are sold.\n    * family identifies the type of product sold.\n    * sales gives the total sales for a product family at a particular store at a given date. Fractional values are possible since products can be sold in   fractional units (1.5 kg of cheese, for instance, as opposed to 1 bag of chips).\n    * onpromotion gives the total number of items in a product family that were being promoted at a store at a given date.\n\n2. **test.csv**\nThe test data, having the same features as the training data. Predict the target sales for the dates in this file. The dates in the test data are for the 15 days after the last date in the training data.\n\n3. **sample_submission.csv**\nA sample submission file in the correct format.\n\n4. **stores.csv**\nStore metadata, including city, state, type, and cluster.\n\ncluster is a grouping of similar stores.\n\n5. **oil.csv**\nDaily oil price. Includes values during both the train and test data timeframes. (Ecuador is an oil-dependent country and it's economical health is highly vulnerable to shocks in oil prices.)\n\n6. **holidays_events.csv**\nHolidays and Events, with metadata\nNOTE: Pay special attention to the transferred column. A holiday that is transferred officially falls on that calendar day, but was moved to another date by the government. A transferred day is more like a normal day than a holiday. To find the day that it was actually celebrated, look for the corresponding row where type is Transfer. For example, the holiday Independencia de Guayaquil was transferred from 2012-10-09 to 2012-10-12, which means it was celebrated on 2012-10-12. Days that are type Bridge are extra days that are added to a holiday (e.g., to extend the break across a long weekend). These are frequently made up by the type Work Day which is a day not normally scheduled for work (e.g., Saturday) that is meant to payback the Bridge.\nAdditional holidays are days added a regular calendar holiday, for example, as typically happens around Christmas (making Christmas Eve a holiday).\n\n7. **Additional Notes**\nWages in the public sector are paid every two weeks on the 15 th and on the last day of the month. Supermarket sales could be affected by this.\nA magnitude 7.8 earthquake struck Ecuador on April 16, 2016. People rallied in relief efforts donating water and other first need products which greatly affected supermarket sales for several weeks after the earthquake.","metadata":{}},{"cell_type":"code","source":"sales = pd.read_csv('../input/store-sales-time-series-forecasting/train.csv',\n                          dtype={'store_nbr': 'category', 'family': 'category', 'sales': 'float32'}, \n                          parse_dates=['date'],infer_datetime_format=True)\nsales.head().style.set_properties(**{'background-color': 'LavenderBlush'})","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:05.125526Z","iopub.execute_input":"2022-07-27T08:03:05.126082Z","iopub.status.idle":"2022-07-27T08:03:07.874615Z","shell.execute_reply.started":"2022-07-27T08:03:05.126009Z","shell.execute_reply":"2022-07-27T08:03:07.873171Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sales.tail().style.set_properties(**{'background-color': 'LavenderBlush'})","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:07.884522Z","iopub.execute_input":"2022-07-27T08:03:07.885395Z","iopub.status.idle":"2022-07-27T08:03:07.906841Z","shell.execute_reply.started":"2022-07-27T08:03:07.885351Z","shell.execute_reply":"2022-07-27T08:03:07.905342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sales.store_nbr.nunique(), sales.family.nunique(),  sales.date.nunique()","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:07.909307Z","iopub.execute_input":"2022-07-27T08:03:07.909757Z","iopub.status.idle":"2022-07-27T08:03:07.977952Z","shell.execute_reply.started":"2022-07-27T08:03:07.909710Z","shell.execute_reply":"2022-07-27T08:03:07.976794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Rearrange the dataset in such a way as to determine the Category sales per store.","metadata":{}},{"cell_type":"code","source":"sales_df = sales.drop('onpromotion', axis = 1).pivot_table(index=['date'], values = ['sales'],columns = ['store_nbr','family'],fill_value = 0)\nsales_df.columns = [\"_\".join(x) for x in sales_df.columns.ravel()]\nsales_df.head().style.set_properties(**{'background-color': 'LavenderBlush'})","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:07.979555Z","iopub.execute_input":"2022-07-27T08:03:07.979919Z","iopub.status.idle":"2022-07-27T08:03:13.590398Z","shell.execute_reply.started":"2022-07-27T08:03:07.979887Z","shell.execute_reply":"2022-07-27T08:03:13.589217Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To confirm whether the dataset contains any missing values and/or duplicates create a complete list of DateTime with 15 minutes interval from the starting to end point and check whether it matches with the index list of dataset.","metadata":{}},{"cell_type":"code","source":"# sort by dates\nsales_df.sort_index(inplace = True)\n\n#creating datetime list with boundaries of raw data series, hourly frequency\ndatelist = pd.date_range(datetime.datetime(2013,1,1), datetime.datetime(2017,8,15), freq='D').tolist()\n\n#extracting raw data series indices\nidx_list = sales_df.index.to_list()\n\n#checking for anomalies by comparing the two\nprint(idx_list == datelist)\n#searching for anomalies\nprint(\"\\n No. of elements in full list:\", len(datelist), \"\\n No. of indices:\", len(idx_list), \"\\n No. of elements in set of indices:\", len(set(idx_list)))","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:13.592037Z","iopub.execute_input":"2022-07-27T08:03:13.592472Z","iopub.status.idle":"2022-07-27T08:03:13.614129Z","shell.execute_reply.started":"2022-07-27T08:03:13.592431Z","shell.execute_reply":"2022-07-27T08:03:13.613146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Number of indices in the dataframe is lesser than that should be in the date range. Notice there are two missing dates in the index. Add the missing date in the DatetimeIndex by replacing it with a new index using reindex(). It sets NaN to values whose row or column label is new.","metadata":{}},{"cell_type":"code","source":"sales_df = sales_df.reindex(pd.date_range(datetime.datetime(2013,1,1), datetime.datetime(2017,8,15), freq='D'), fill_value=0)","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:13.615480Z","iopub.execute_input":"2022-07-27T08:03:13.616010Z","iopub.status.idle":"2022-07-27T08:03:13.634423Z","shell.execute_reply.started":"2022-07-27T08:03:13.615977Z","shell.execute_reply":"2022-07-27T08:03:13.633340Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Check again to make sure all dates are present in the dataframe. ","metadata":{}},{"cell_type":"code","source":"# sort by dates\nsales_df.sort_index(inplace = True)\n\n#creating datetime list with boundaries of raw data series, hourly frequency\ndatelist = pd.date_range(datetime.datetime(2013,1,1), datetime.datetime(2017,8,15), freq='D').tolist()\n\n#extracting raw data series indices\nidx_list = sales_df.index.to_list()\n\n#checking for anomalies by comparing the two\nprint(idx_list == datelist)\n#searching for anomalies\nprint(\"\\n No. of elements in full list:\", len(datelist), \"\\n No. of indices:\", len(idx_list), \"\\n No. of elements in set of indices:\", len(set(idx_list)))","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:13.635606Z","iopub.execute_input":"2022-07-27T08:03:13.636386Z","iopub.status.idle":"2022-07-27T08:03:13.657194Z","shell.execute_reply.started":"2022-07-27T08:03:13.636351Z","shell.execute_reply":"2022-07-27T08:03:13.656367Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **<span style = 'color:green'>Visualize the Time Series to study the original plot</span>**<a id ='Visualize'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\n\n### **<span style = 'color:brown'>Sale Categories of Store Number 1</span>**<a id = 'Store1'></a>","metadata":{}},{"cell_type":"markdown","source":"sales_df2 = sales.pivot(index ='date', values = ['sales'],columns = ['store_nbr','family'])\nsales_df2.style.set_properties(**{'background-color': 'LavenderBlush'})","metadata":{"execution":{"iopub.status.busy":"2022-07-26T20:54:55.739094Z","iopub.execute_input":"2022-07-26T20:54:55.739515Z","iopub.status.idle":"2022-07-26T20:58:09.847035Z","shell.execute_reply.started":"2022-07-26T20:54:55.739483Z","shell.execute_reply":"2022-07-26T20:58:09.845526Z"}}},{"cell_type":"code","source":"import plotly.express as px\nfig = px.line(sales_df, x=sales_df.index, y=sales_df.columns[0:32],\n              title='Sale Categories of Store Number 1', width=1500)\n\nfig.update_layout(\n    updatemenus=[\n        dict(\n            active=0,\n            buttons=list([dict(label=\"All\",\n                     method=\"update\",\n                     args=[{\"visible\": [True for _ in range(33)]},\n                           {\"title\": \"Sale Categories of Store Number 1\",\n                            \"annotations\": []}])]) + list([\n                dict(label=f\"{j}\",\n                     method=\"update\",\n                     args=[{\"visible\": [True if i==idx else False for i in range(33)]},\n                           {\"title\": f\"{j}\",\n                            \"annotations\": []}]) for idx,j in enumerate(sales_df.columns[0:32])])\n            )])\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:13.658566Z","iopub.execute_input":"2022-07-27T08:03:13.659099Z","iopub.status.idle":"2022-07-27T08:03:15.406337Z","shell.execute_reply.started":"2022-07-27T08:03:13.659042Z","shell.execute_reply":"2022-07-27T08:03:15.404929Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Total Sales Per Year</span>**<a id = 'Total'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)","metadata":{}},{"cell_type":"code","source":"Total_sales = sales.groupby('date').sum().squeeze()\nTotal_sales = Total_sales.reindex(pd.date_range(datetime.datetime(2013,1,1), datetime.datetime(2017,8,15), freq='D'), fill_value=0)\nTotal_sales.head().style.set_properties(**{'background-color': 'LavenderBlush'})","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:15.407733Z","iopub.execute_input":"2022-07-27T08:03:15.408660Z","iopub.status.idle":"2022-07-27T08:03:15.566599Z","shell.execute_reply.started":"2022-07-27T08:03:15.408619Z","shell.execute_reply":"2022-07-27T08:03:15.565791Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.express as px\nfig = px.line(Total_sales, x=Total_sales.index, y=Total_sales.columns[1],\n              title='Total Sales per year', width=1300)\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:15.568333Z","iopub.execute_input":"2022-07-27T08:03:15.569027Z","iopub.status.idle":"2022-07-27T08:03:15.685413Z","shell.execute_reply.started":"2022-07-27T08:03:15.568980Z","shell.execute_reply":"2022-07-27T08:03:15.684268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Daily Store Sales Data</span>**<a id = 'Daily'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)","metadata":{}},{"cell_type":"code","source":"sales = pd.read_csv('../input/store-sales-time-series-forecasting/train.csv', \n                    dtype={'store_nbr': 'int32', 'family': 'category', 'sales': 'float32'}, \n                          parse_dates=['date'],infer_datetime_format=True)\nstore_sales=sales.sort_values(by ='store_nbr')\nstore_sales['store_nbr'] = store_sales['store_nbr'].astype('str')\nstore_sales = store_sales.groupby(['date', 'store_nbr']).mean()\nstore_sales = store_sales.drop('onpromotion', axis = 1).pivot_table(index=['date'], values = ['sales'],columns = ['store_nbr'],fill_value = 0)\nstore_sales = store_sales.reindex(pd.date_range(datetime.datetime(2013,1,1), datetime.datetime(2017,8,15), freq='D'), fill_value=0)\nstore_sales.columns = [\"_\".join(x) for x in store_sales.columns.ravel()]\nstore_sales","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:15.687037Z","iopub.execute_input":"2022-07-27T08:03:15.687429Z","iopub.status.idle":"2022-07-27T08:03:22.543443Z","shell.execute_reply.started":"2022-07-27T08:03:15.687394Z","shell.execute_reply":"2022-07-27T08:03:22.541929Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.express as px\nfig = px.line(store_sales, x=store_sales.index, y=store_sales.columns[0:],\n              title='Daily Store Sales Data', width=1300)\n\nfig.update_layout(\n    updatemenus=[\n        dict(\n            active=0,\n            buttons=list([dict(label=\"All\",\n                     method=\"update\",\n                     args=[{\"visible\": [True for _ in range(54)]},\n                           {\"title\": \"Daily Store Sales Data\",\n                            \"annotations\": []}])]) + list([\n                dict(label=f\"{j}\",\n                     method=\"update\",\n                     args=[{\"visible\": [True if i==idx else False for i in range(54)]},\n                           {\"title\": f\"{j}\",\n                            \"annotations\": []}]) for idx,j in enumerate(store_sales.columns[0:])])\n            )])\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:22.545050Z","iopub.execute_input":"2022-07-27T08:03:22.545452Z","iopub.status.idle":"2022-07-27T08:03:25.573518Z","shell.execute_reply.started":"2022-07-27T08:03:22.545419Z","shell.execute_reply":"2022-07-27T08:03:25.571476Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **<span style = 'color:green'>Feature Engineering & Data Preprocessing </span>**<a id ='Preprocessing'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\nCode Credits: [Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}},{"cell_type":"code","source":"class DataSet:\n    \"\"\"\n    Preprocessing.\n    \"\"\"\n    def __init__(self, horizon, back_horizon):\n        self.horizon = horizon\n        self.back_horizon = back_horizon\n    \n    def preprocessing(self, y, date, train_size=0.7, val_size=0.2):\n        \n        y = y.copy().astype('float')\n\n        train = y[:int(train_size*len(y))]\n        val = y[int(train_size*len(y))-self.back_horizon:int((train_size+val_size)*len(y))]\n        test = y[int((train_size+val_size)*len(y))-self.back_horizon:]\n        train_date = date[:int(train_size*len(y))]\n        val_date = date[int(train_size*len(y))-self.back_horizon:int((train_size+val_size)*len(y))]\n        test_date = date[int((train_size+val_size)*len(y))-self.back_horizon:]\n\n        # Training set\n        self.X_train, self.y_train, self.train_date = self.create_sequences(train, \n                                                                            train, \n                                                                            train_date,\n                                                                            self.horizon, \n                                                                            self.back_horizon)\n        # Validation set\n        self.X_val, self.y_val, self.val_date = self.create_sequences(val,\n                                                                      val,\n                                                                      val_date,\n                                                                      self.horizon,\n                                                                      self.back_horizon)\n        # Testing set\n        self.X_test, self.y_test, self.test_date = self.create_sequences(test,\n                                                                         test,\n                                                                         test_date,\n                                                                         self.horizon,\n                                                                         self.back_horizon)\n\n        # training on all database\n        self.X_train_all, self.y_train_all, self.train_all_date = self.create_sequences(y, \n                                                                                        y,\n                                                                                        date,\n                                                                                        self.horizon,\n                                                                                        self.back_horizon)\n            \n    @staticmethod\n    def create_sequences(X, y, d, horizon, time_steps):\n        Xs, ys, ds = [], [], []\n        for col in range(X.shape[1]):\n            for i in range(0, len(X)-time_steps-horizon, 1):\n                Xs.append(X[i:(i+time_steps), col])\n                ys.append(y[(i+time_steps):(i+time_steps+horizon), col])\n                ds.append(d[(i+time_steps):(i+time_steps+horizon)])\n\n        return np.array(Xs), np.array(ys), np.array(ds)","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:25.576286Z","iopub.execute_input":"2022-07-27T08:03:25.576981Z","iopub.status.idle":"2022-07-27T08:03:25.601700Z","shell.execute_reply.started":"2022-07-27T08:03:25.576910Z","shell.execute_reply":"2022-07-27T08:03:25.600127Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **<span style='color:green'>Define Interpretable NBeats model</span>**<a id ='Build'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\nNeural Basis Expansion Analysis Time Series (NBEATS) is based on backward and forward residual links and a very deep stack of fully-connected layers. Model type can be ‘generic’ or ‘interpretable’. Interpretable models which describes series as a function of trend and seasonality, can help to better understand interactions between variables also.\n\n### **<span style = 'color:brown'>Code to structure NBeats architecture</span>**<a id ='Architecture'></a>\n\n\nCode Credits: [Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}},{"cell_type":"code","source":"import tensorflow as tf\nclass Stack(tf.keras.layers.Layer):\n\n    \"\"\"A stack is a series of blocks where each block produce two outputs, the horizon and the back_horizon. \n    All of the outputs are sum up which compose the stack output while each residual back_horizon is given to the following block.\n    \n    Parameters\n    ----------\n    blocks: list of `TrendBlock`, `SeasonalityBlock` or `GenericBlock`.\n        Define blocks in a stack.\n    \"\"\"\n    def __init__(self, blocks, **kwargs):\n        \n        super().__init__(**kwargs)\n\n        self._blocks = blocks\n\n    def call(self, inputs):\n\n        y_horizon = 0.\n        for block in self._blocks:\n            residual_y, y_back_horizon = block(inputs) # shape: (n_quantiles, Batch_size, horizon), (Batch_size, back_horizon)\n            inputs = tf.subtract(inputs, y_back_horizon)\n            y_horizon = tf.add(y_horizon, residual_y) # shape: (n_quantiles, Batch_size, horizon)\n\n        return y_horizon, inputs","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:25.603351Z","iopub.execute_input":"2022-07-27T08:03:25.603767Z","iopub.status.idle":"2022-07-27T08:03:25.618506Z","shell.execute_reply.started":"2022-07-27T08:03:25.603729Z","shell.execute_reply":"2022-07-27T08:03:25.617350Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Define NBEATS wrapper to create generic model and interpretable model</span>**<a id ='NBEATS-wrapper'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\n\nCode Credits: [Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}},{"cell_type":"code","source":"class N_BEATS(tf.keras.Model):\n    \"\"\"This class compute the N-BEATS model. This is a univariate model which can be\n     interpretable or generic. It's strong advantage is its internal structure which allows us \n     to extract the trend and the seasonality of a temporal serie. It's available from the attributes\n     `seasonality` and `trend`. This is an unofficial implementation.\n\n     `@inproceedings{\n        Oreshkin2020:N-BEATS,\n        title={{N-BEATS}: Neural basis expansion analysis for interpretable time series horizoning},\n        author={Boris N. Oreshkin and Dmitri Carpov and Nicolas Chapados and Yoshua Bengio},\n        booktitle={International Conference on Learning Representations},\n        year={2020},\n        url={https://openreview.net/forum?id=r1ecqn4YwB}\n        }`\n    \n    Parameter\n    ---------\n    stacks: list of `Stack` layer.\n        Define the stack to use in nbeats model. It can be full of `TrendBlock`, `SeasonalityBlock` or `GenereicBlock`.\n    \"\"\"\n    def __init__(self, \n                 stacks,\n                 **kwargs):\n                \n        super().__init__(**kwargs)\n\n        self._stacks = stacks\n\n    def call(self, inputs):\n        self._residuals_y = tf.TensorArray(tf.float32, size=len(self._stacks)) # Stock trend and seasonality curves during inference\n        y_horizon = 0.\n        for idx, stack in enumerate(self._stacks):\n            residual_y, inputs = stack(inputs)\n            self._residuals_y.write(idx, residual_y)\n            y_horizon = tf.add(y_horizon, residual_y)\n\n        return y_horizon\n\n    @property\n    def seasonality(self):\n        return self._residuals_y.stack()[1]\n\n    @property\n    def trend(self):\n        return self._residuals_y.stack()[0]","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:25.620298Z","iopub.execute_input":"2022-07-27T08:03:25.620965Z","iopub.status.idle":"2022-07-27T08:03:25.634451Z","shell.execute_reply.started":"2022-07-27T08:03:25.620930Z","shell.execute_reply":"2022-07-27T08:03:25.632658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Define Generic Block</span>**<a id ='Generic'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\n\nCode Credits: [Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}},{"cell_type":"code","source":"from tensorflow import keras\nclass GenericBlock(keras.layers.Layer):\n    \"\"\"\n    Generic block definition as described in the paper. \n    We can't have explanation from this kind of block because g coefficients are learnt.\n    \n    Parameter\n    ---------\n    horizon: integer\n        Horizon time to horizon.\n    back_horizon: integer\n        Past to rebuild.\n    nb_neurons: integer\n        Number of neurons in Fully connected layers.\n    back_neurons: integer\n        Number of back_horizon expansion coefficients.\n    fore_neurons: integer\n        Number of horizon expansion coefficients.\n    \"\"\"\n    def __init__(self, horizon, \n                       back_horizon, \n                       n_neurons, \n                       n_quantiles,\n                       dropout_rate,\n                        **kwargs):\n        \n        super().__init__(**kwargs)\n        \n        self._FC_stack = [keras.layers.Dense(n_neurons, \n                                            activation='relu', \n                                            kernel_initializer=\"glorot_uniform\") for _ in range(4)]\n        \n        self._dropout = tf.keras.layers.Dropout(dropout_rate)\n\n        \n        self._FC_back_horizon = self.add_weight(shape=(n_neurons, n_neurons), \n                                           trainable=True,\n                                           initializer=\"glorot_uniform\",\n                                           name='FC_back_horizon_generic')\n\n        self._FC_horizon = self.add_weight(shape=(n_quantiles, n_neurons, n_neurons), \n                                           trainable=True,\n                                           initializer=\"glorot_uniform\",\n                                           name='FC_horizon_generic')\n        \n        self._back_horizon = keras.layers.Dense(back_horizon, \n                                           kernel_initializer=\"glorot_uniform\")\n        \n        self._horizon = keras.layers.Dense(horizon, \n                                           kernel_initializer=\"glorot_uniform\")\n        \n    def call(self, inputs):\n        # shape: (Batch_size, back_horizon)\n        for dense_layer in self._FC_stack:\n            inputs = dense_layer(inputs) # shape: (Batch_size, nb_neurons)\n            inputs = self._dropout(inputs, training=True) # We bind first layers by a dropout \n            \n        theta_horizon = inputs @ self._FC_horizon # shape: (n_quantiles, Batch_size, 2 * fourier order)\n        theta_back_horizon = inputs @ self._FC_back_horizon # shape: (Batch_size, 2 * fourier order)\n        \n        y_back_horizon = self._back_horizon(theta_back_horizon) # shape: (Batch_size, back_horizon)\n        y_horizon = self._horizon(theta_horizon) # shape: (n_quantiles, Batch_size, horizon)\n        \n        return y_horizon, y_back_horizon","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:25.636356Z","iopub.execute_input":"2022-07-27T08:03:25.637345Z","iopub.status.idle":"2022-07-27T08:03:25.657438Z","shell.execute_reply.started":"2022-07-27T08:03:25.637292Z","shell.execute_reply":"2022-07-27T08:03:25.656134Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Define Trend Block</span>**<a id ='Trend'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\n\nCode Credits: [Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}},{"cell_type":"code","source":"class TrendBlock(tf.keras.layers.Layer):\n    \"\"\" Trend block definition. Output layers are constrained which define polynomial function of small degree p.\n    Therefore it is possible to get explanation from this block.\n    \n    Parameter\n    ---------\n    p_degree: integer\n        Degree of the polynomial function.\n    horizon: integer\n        Horizon time to horizon.\n    back_horizon: integer\n        Past to rebuild.\n    n_neurons: integer\n        Number of neurons in Fully connected layers.\n    n_quantiles: Integer.\n        Number of quantiles in `QuantileLossError`.\n    \"\"\"\n    def __init__(self, \n                 horizon, \n                 back_horizon,\n                 p_degree,   \n                 n_neurons, \n                 n_quantiles, \n                 dropout_rate,\n                 **kwargs):\n\n        super().__init__(**kwargs)\n        \n        self._p_degree = tf.reshape(tf.range(p_degree + 1, dtype='float32'), shape=(-1, 1)) # Shape (-1, 1) in order to broadcast horizon to all p degrees\n        self._horizon = tf.cast(horizon, dtype='float32') \n        self._back_horizon = tf.cast(back_horizon, dtype='float32')\n        self._n_neurons = n_neurons \n        self._n_quantiles = n_quantiles\n\n        self._FC_stack = [tf.keras.layers.Dense(n_neurons, \n                                            activation='relu', \n                                            kernel_initializer=\"glorot_uniform\") for _ in range(4)]\n        \n        self._dropout = tf.keras.layers.Dropout(dropout_rate)\n        \n        self._FC_back_horizon = self.add_weight(shape=(n_neurons, p_degree + 1), \n                                           trainable=True,\n                                           initializer=\"glorot_uniform\",\n                                           name='FC_back_horizon_trend')\n        \n        self._FC_horizon = self.add_weight(shape=(n_quantiles, n_neurons, p_degree + 1),\n                                           trainable=True,\n                                           initializer=\"glorot_uniform\",\n                                           name='FC_horizon_trend')\n\n        self._horizon_coef = (tf.range(self._horizon) / self._horizon) ** self._p_degree\n        self._back_horizon_coef = (tf.range(self._back_horizon) / self._back_horizon) ** self._p_degree\n        \n    def call(self, inputs):\n\n        for dense in self._FC_stack:\n            inputs = dense(inputs) # shape: (Batch_size, n_neurons)\n            inputs = self._dropout(inputs, training=True) # We bind first layers by a dropout \n            \n        theta_back_horizon = inputs @ self._FC_back_horizon # shape: (Batch_size, p_degree)\n        theta_horizon = inputs @ self._FC_horizon # shape: (n_quantiles, Batch_size, p_degree)\n\n        y_back_horizon = theta_back_horizon @ self._back_horizon_coef # shape: (Batch_size, back_horizon)\n        y_horizon = theta_horizon @ self._horizon_coef # shape: (n_quantiles, Batch_size, horizon)\n        \n        return y_horizon, y_back_horizon\n    ","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:25.659703Z","iopub.execute_input":"2022-07-27T08:03:25.660666Z","iopub.status.idle":"2022-07-27T08:03:25.679201Z","shell.execute_reply.started":"2022-07-27T08:03:25.660626Z","shell.execute_reply":"2022-07-27T08:03:25.677672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Define Seasonality Block</span>**<a id ='Seasonality'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\n\nCode Credits: [Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}},{"cell_type":"code","source":"class SeasonalityBlock(tf.keras.layers.Layer):\n    \"\"\"Seasonality block definition. Output layers are constrained which define fourier series. \n    Each expansion coefficent then become a coefficient of the fourier serie. As each block and each \n    stack outputs are sum up, we decided to introduce fourier order and multiple seasonality periods.\n    Therefore it is possible to get explanation from this block.\n    \n    Parameters\n    ----------\n    horizon: integer\n        Horizon time to horizon.\n    back_horizon: integer\n        Past to rebuild.\n    n_neurons: integer\n        Number of neurons in Fully connected layers.\n    periods: Integer.\n        fourier serie period. The paper set this parameter to `horizon/2`.\n    back_periods: Integer.\n        fourier serie back period. The paper set this parameter to `back_horizon/2`.\n    horizon_fourier_order: Integer.\n        Higher values signifies complex fourier serie\n    back_horizon_fourier_order: Integer.\n        Higher values signifies complex fourier serie\n    n_quantiles: Integer.\n        Number of quantiles in `QuantileLossError`.\n    \"\"\"\n    def __init__(self,\n                 horizon,\n                 back_horizon,\n                 n_neurons, \n                 periods, \n                 back_periods, \n                 horizon_fourier_order,\n                 back_horizon_fourier_order,\n                 n_quantiles,\n                 dropout_rate,\n                 **kwargs):\n        \n        super().__init__(**kwargs)\n\n        self._horizon = horizon\n        self._back_horizon = back_horizon\n        self._periods = tf.cast(tf.reshape(periods, (1, -1)), 'float32') # Broadcast horizon on multiple periods\n        self._back_periods = tf.cast(tf.reshape(back_periods, (1, -1)), 'float32')  # Broadcast back horizon on multiple periods\n        self._horizon_fourier_order = tf.reshape(tf.range(horizon_fourier_order, dtype='float32'), shape=(-1, 1)) # Broadcast horizon on multiple fourier order\n        self._back_horizon_fourier_order = tf.reshape(tf.range(back_horizon_fourier_order, dtype='float32'), shape=(-1, 1)) # Broadcast horizon on multiple fourier order\n\n        # Workout the number of neurons needed to compute seasonality coefficients\n        horizon_neurons = tf.reduce_sum(2 * horizon_fourier_order)\n        back_horizon_neurons = tf.reduce_sum(2 * back_horizon_fourier_order)\n        \n        self._FC_stack = [tf.keras.layers.Dense(n_neurons, \n                                               activation='relu', \n                                               kernel_initializer=\"glorot_uniform\") for _ in range(4)]\n        \n        self._dropout = tf.keras.layers.Dropout(dropout_rate)   \n        \n        self._FC_back_horizon = self.add_weight(shape=(n_neurons, back_horizon_neurons), \n                                           trainable=True,\n                                           initializer=\"glorot_uniform\",\n                                           name='FC_back_horizon_seasonality')\n\n        self._FC_horizon = self.add_weight(shape=(n_quantiles, n_neurons, horizon_neurons), \n                                           trainable=True,\n                                           initializer=\"glorot_uniform\",\n                                           name='FC_horizon_seasonality')\n        \n        # Workout cos and sin seasonality coefficents\n        time_horizon = tf.range(self._horizon, dtype='float32') / self._periods\n        horizon_seasonality = 2 * np.pi * self._horizon_fourier_order * time_horizon\n        self._horizon_coef = tf.concat((tf.cos(horizon_seasonality), \n                                          tf.sin(horizon_seasonality)), axis=0)\n\n        time_back_horizon = tf.range(self._back_horizon, dtype='float32') / self._back_periods\n        back_horizon_seasonality = 2 * np.pi * self._back_horizon_fourier_order * time_back_horizon\n        self._back_horizon_coef = tf.concat((tf.cos(back_horizon_seasonality), \n                                        tf.sin(back_horizon_seasonality)), axis=0)\n        \n    def call(self, inputs):\n\n        for dense in self._FC_stack:\n            inputs = dense(inputs) # shape: (Batch_size, nb_neurons)\n            inputs = self._dropout(inputs, training=True) # We bind first layers by a dropout \n\n        theta_horizon = inputs @ self._FC_horizon # shape: (n_quantiles, Batch_size, 2 * fourier order)\n        theta_back_horizon = inputs @ self._FC_back_horizon # shape: (Batch_size, 2 * fourier order)\n        \n        y_horizon = theta_horizon @ self._horizon_coef # shape: (n_quantiles, Batch_size, horizon)\n        y_back_horizon = theta_back_horizon @ self._back_horizon_coef # shape: (Batch_size, back_horizon)\n    \n        return y_horizon, y_back_horizon","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:25.681288Z","iopub.execute_input":"2022-07-27T08:03:25.682438Z","iopub.status.idle":"2022-07-27T08:03:25.711609Z","shell.execute_reply.started":"2022-07-27T08:03:25.682389Z","shell.execute_reply":"2022-07-27T08:03:25.710659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Define Interpretable Block</span>**<a id ='Interpretable'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\nCode Credits: [Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}},{"cell_type":"code","source":"def create_interpretable_nbeats(horizon, \n                                   back_horizon,\n                                   p_degree,   \n                                   trend_n_neurons, \n                                   seasonality_n_neurons, \n                                   periods, \n                                   back_periods, \n                                   horizon_fourier_order,\n                                   back_horizon_fourier_order,\n                                   n_quantiles, \n                                   share=True,\n                                   dropout_rate=0.1,\n                                   **kwargs):\n    \n    \"\"\"Wrapper to create interpretable model. check nbeats doc to know more about parameters.\"\"\"\n    \n    if share is True:\n        trend_block = TrendBlock(horizon=horizon, \n                                 back_horizon=back_horizon, \n                                 p_degree=p_degree, \n                                 n_neurons=trend_n_neurons, \n                                 n_quantiles=n_quantiles,\n                                 dropout_rate=dropout_rate,\n                                 **kwargs)\n        \n        seasonality_block = SeasonalityBlock(horizon=horizon, \n                                 back_horizon=back_horizon, \n                                 periods=periods, \n                                 back_periods=back_periods, \n                                 horizon_fourier_order=horizon_fourier_order,\n                                 back_horizon_fourier_order = back_horizon_fourier_order,\n                                 n_neurons=seasonality_n_neurons, \n                                 n_quantiles=n_quantiles,dropout_rate=dropout_rate, **kwargs)\n        \n        trendblocks = [trend_block for _ in range(3)]\n        seasonalityblocks = [seasonality_block for _ in range(3)]\n    else:\n        trendblocks = [TrendBlock(horizon=horizon, \n                                 back_horizon=back_horizon, \n                                 p_degree=p_degree, \n                                 n_neurons=trend_n_neurons, \n                                 n_quantiles=n_quantiles, \n                                  dropout_rate=dropout_rate,\n                                 **kwargs) for _ in range(3)]\n        seasonalityblocks = [SeasonalityBlock(horizon=horizon, \n                                 back_horizon=back_horizon, \n                                 periods=periods, \n                                 back_periods=back_periods, \n                                 horizon_fourier_order=horizon_fourier_order,\n                                 back_horizon_fourier_order = back_horizon_fourier_order,\n                                 n_neurons=seasonality_n_neurons, \n                                 n_quantiles=n_quantiles,dropout_rate=dropout_rate, **kwargs) for _ in range(3)]\n        \n    trendstacks = Stack(trendblocks)\n    seasonalitystacks = Stack(seasonalityblocks)\n    \n    return N_BEATS([trendstacks, seasonalitystacks])\n\ndef create_generic_nbeats(horizon,\n                          back_horizon, \n                           n_neurons, \n                           n_quantiles,\n                           n_blocks,\n                           n_stacks,\n                           share=True,\n                           dropout_rate=0.1,\n                           **kwargs):\n\n    \"\"\"Wrapper to create interpretable model. check nbeats doc to know more about parameters.\"\"\"\n    generic_stacks = []\n    if share is True:\n        for stack in range(n_stacks):\n            generic_block = GenericBlock(horizon=horizon, \n                              back_horizon=back_horizon, \n                               n_neurons=n_neurons, \n                               n_quantiles=n_quantiles,\n                               dropout_rate=0.1,\n                               **kwargs)\n\n            generic_blocks = [generic_block for _ in range(n_blocks)]\n            generic_stacks.append(Stack(generic_blocks))\n            \n    else:\n        for stack in range(n_stacks):\n            generic_blocks = [GenericBlock(horizon=horizon, \n                              back_horizon=back_horizon, \n                               n_neurons=n_neurons, \n                               n_quantiles=n_quantiles, \n                               dropout_rate=0.1,\n                               **kwargs) for _ in range(n_blocks)]\n            \n            generic_stacks.append(Stack(generic_blocks))\n    \n    return N_BEATS(generic_stacks)","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:25.712993Z","iopub.execute_input":"2022-07-27T08:03:25.713550Z","iopub.status.idle":"2022-07-27T08:03:25.733601Z","shell.execute_reply.started":"2022-07-27T08:03:25.713517Z","shell.execute_reply":"2022-07-27T08:03:25.732035Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Define quantile loss error function</span>**<a id ='Loss'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\nCode Credits: [Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}},{"cell_type":"code","source":"\"\"\"Functions to assess loss.\"\"\"\nfrom tensorflow.python.keras.losses import LossFunctionWrapper\nfrom tensorflow.python.keras.utils import losses_utils\ndef quantile_loss(y_true, y_pred, quantiles):\n    \n    \"\"\"Calculate the quantile loss function, summed across all quantile outputs.\n\n    Parameters\n    ----------\n    y_true : ndarray or dataframe or list or Tensor of shape `[batch_size, d0, .. dN]`.\n        Ground truth values.\n\n    y_pred : ndarray or dataframe or list or Tensor of shape `[batch_size, d0, .. dN]`.\n        The predicted values.\n\n    n_quantiles : ndarray or dataframe or list or Tensor of shape `[batch_size, d0, .. dN]`.\n        The set of output n_quantiles on which is calculated the quantile loss.\n\n    Returns\n    -------\n    error : Tensor of shape `[batch_size, d0, .. dN]`.\n        The error in %.\n    \"\"\"\n\n    y_true = tf.cast(y_true, dtype=tf.float32)\n    y_pred = tf.cast(y_pred, dtype=tf.float32)\n    quantiles = tf.convert_to_tensor(quantiles)\n    diff = tf.transpose(y_true - y_pred)\n\n    quantile_loss = (quantiles * tf.clip_by_value(diff, 0., np.inf) +\n                    (1 - quantiles) * tf.clip_by_value(-diff, 0., np.inf))\n    \n    M = tf.cast(tf.shape(y_true)[0], dtype=tf.float32)\n    error = quantile_loss / M\n    \n    sum_quantiles = tf.reduce_sum(error, axis=-1)\n    return tf.reduce_sum(sum_quantiles, axis=tf.range(sum_quantiles.shape.rank-1))\n\nclass QuantileLossError(LossFunctionWrapper):\n    \"\"\"Calculate the quantile loss error between `y_true`and `y_pred` across all examples.\n    Standalone usage:\n    >>> y_true = [[0., 1., 2.], [0., 0., 4.]]\n    >>> y_pred = [[1., 1., 2.], [1., 0., 3.]]\n    >>> # Using 'auto'/'sum_over_batch_size' reduction type.\n    >>> ql = QuantileLossError(quantiles=[0.5])\n    >>> ql(y_true, y_pred).numpy()\n    0.5\n    >>> # Calling with 'sample_weight'.\n    >>> ql(y_true, y_pred, sample_weight=[0.7, 0.3]).numpy()\n    0.25\n    >>> # Using 'AUTO' reduction type.\n    >>> ql = QuantileLossError(quantiles=[0.5],\n    ...     reduction=tf.keras.losses.Reduction.AUTO)\n    >>> ql(y_true, y_pred).numpy()\n    0.25\n    >>> # Using 'none' reduction type.\n    >>> ql = QuantileLossError(quantiles=[0.5],\n    ...     reduction=tf.keras.losses.Reduction.NONE)\n    >>> ql(y_true, y_pred).numpy()\n    array([0.25, 0.25], dtype=float32)\n    >>> # Using multiple quantiles.\n    >>> ql = QuantileLossError(quantiles=[0.1, 0.5, 0.9])\n    >>> ql(y_true, y_pred).numpy()\n    1.5\n\n    Usage with the `compile()` API:\n    ```python\n    model.compile(optimizer='sgd', loss=tf.keras.losses.QuantileLossError())\n    ```\n    \"\"\"\n\n    def __init__(self,\n                 quantiles,\n                 reduction=losses_utils.ReductionV2.SUM,\n                 name='quantile_loss'):\n        super(QuantileLossError, self).__init__(\n            quantile_loss, quantiles=quantiles, name=name, reduction=reduction)","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:25.735166Z","iopub.execute_input":"2022-07-27T08:03:25.735873Z","iopub.status.idle":"2022-07-27T08:03:25.753449Z","shell.execute_reply.started":"2022-07-27T08:03:25.735838Z","shell.execute_reply":"2022-07-27T08:03:25.752139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Define plot implementation function</span>**<a id ='Plot'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\nCode Credits: [Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}},{"cell_type":"code","source":"def confidence_interval(time_series, quantile):\n    \n    mean = tf.reduce_mean(time_series, axis=0)\n    standard_deviation = tf.math.reduce_std(time_series, axis=0)\n    \n    # (-1, 1, 1, 1) Broadcast std with shape (quantile, timesteps, horizon)\n    max_interval = mean + tf.reshape(tfd.Normal(loc=0, scale=1).quantile(quantile), (-1, 1, 1, 1)) * standard_deviation \n    min_interval = mean - tf.reshape(tfd.Normal(loc=0, scale=1).quantile(quantile), (-1, 1, 1, 1)) * standard_deviation \n    \n    return mean, min_interval, max_interval\n\ndef fig_add_trace(fig, y, x, row, col):\n    for data, fill, name, color, showlegend in y:\n        fig.add_trace(\n            go.Scatter(\n                name=name,\n                x=x,\n                y=data,\n                fill=fill,\n                line=dict(color=color),\n                fillcolor=color,\n                showlegend=showlegend),\n            row=row, col=col)\n        \ndef plot_results_nbeats(y_pred,\n                        date_outputs,\n                        y_true=None,\n                        date_history=None,\n                        seasonality=None,\n                        trend=None):\n    \n    date = pd.to_datetime(date_outputs[0])\n    date_history = date_history if date_history is not None else pd.to_datetime(date_outputs[0])\n    y_pred_mean, y_pred_min_interval, y_pred_max_interval = confidence_interval(y_pred, [0.99])\n    \n    if trend is None:\n        fig = make_subplots(\n        subplot_titles=['True Vs Predicted','Trend','Seasonality', 'Overall trend', 'Overall seasonality'],\n        rows=1, cols=1,\n        vertical_spacing=0.1,\n        horizontal_spacing=0.05,\n        specs=[[{\"type\": \"scatter\"}]])\n        \n        if y_true is not None:\n    \n            # Trace ground truth\n            fig.add_trace(\n                    go.Scatter(\n                        name=\"y_true\",\n                        x=date_history,\n                        y=y_true[0],\n                        line=dict(color=\"green\")),\n                    row=1, col=1\n                )  \n    \n    else:\n        fig = make_subplots(\n        subplot_titles=['True Vs Predicted','Trend','Seasonality', 'Overall trend', 'Overall seasonality'],\n        rows=3, cols=2,\n        vertical_spacing=0.1,\n        horizontal_spacing=0.05,\n        column_widths=[0.8, 0.6],\n        row_heights=[0.8, 0.8, 0.8],\n        specs=[[{\"type\": \"scatter\", \"rowspan\": 2}, {\"type\": \"scatter\"}],\n               [        None      , {\"type\": \"scatter\"}], \n               [{\"type\": \"scatter\", \"colspan\": 2}, None]])\n        \n        if y_true is not None:\n    \n            # Trace ground truth\n            fig.add_trace(\n                    go.Scatter(\n                        name=\"y_true\",\n                        x=date_history,\n                        y=y_true[0],\n                        line=dict(color=\"green\")),\n                    row=1, col=1\n                )  \n            \n        trend_mean, trend_min_interval, trend_max_interval = confidence_interval(trend, [0.99])\n\n        fig_add_trace(fig, zip([trend_min_interval[0, 1, 0], trend_max_interval[0, 1, 0], \n                                trend_mean[1, 0]], \n                              ['none','tonextx', 'none'],\n                              ['Confidence interval', 'Confidence interval', 'Average'],\n                              ['rgba(153,50,204, 0.7)', 'rgba(153,50,204, 0.7)', \n                               'darkorange'], \n                              [ False, False, False]), date, row=1, col=2)\n        \n            # Trace mean\n        fig_add_trace(fig, zip([tf.reduce_mean(trend_min_interval, axis=2)[0, 1], tf.reduce_mean(trend_max_interval, axis=2)[0, 1], \n                            tf.reduce_mean(trend_mean, axis=1)[1]], \n                          ['none','tonextx', 'none'],\n                          ['Confidence interval', 'Confidence interval', 'Average'],\n                          ['rgba(153,50,204, 0.7)', 'rgba(153,50,204, 0.7)', \n                           'darkorange'], \n                          [False, False, False]), date, row=3, col=1)\n    \n    # Trace mean\n    fig_add_trace(fig, zip([y_pred_min_interval[0, 0, 0], y_pred_max_interval[0, 2, 0], \n                            y_pred_min_interval[0, 1, 0], y_pred_max_interval[0, 1, 0], \n                            y_pred_mean[1, 0]], \n                          ['none', 'tonextx', 'none','tonextx', 'none'],\n                          ['Prediction interval', 'Prediction interval', 'Confidence interval', 'Confidence interval', 'Average'],\n                          ['rgba(255, 0, 0, 0.2)', 'rgba(255, 0, 0, 0.2)', \n                           'rgba(153,50,204, 0.7)', 'rgba(153,50,204, 0.7)', \n                           'darkorange'], \n                          [False, True, False, True, True]), date, row=1, col=1)\n    \n    # Trace seasonality\n    if seasonality is not None:\n        seasonality_mean, seasonality_min_interval, seasonality_max_interval = confidence_interval(seasonality, [0.99])\n        fig_add_trace(fig, zip([seasonality_min_interval[0, 1, 0], seasonality_max_interval[0, 1, 0], \n                                seasonality_mean[1, 0]], \n                              ['none','tonextx', 'none'],\n                              ['Confidence interval', 'Confidence interval', 'Average'],\n                              ['rgba(153,50,204, 0.7)', 'rgba(153,50,204, 0.7)', \n                               'darkorange'], \n                              [False, False, False]), date, row=2, col=2)\n    \n    return fig","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:25.755391Z","iopub.execute_input":"2022-07-27T08:03:25.756028Z","iopub.status.idle":"2022-07-27T08:03:25.793258Z","shell.execute_reply.started":"2022-07-27T08:03:25.755883Z","shell.execute_reply":"2022-07-27T08:03:25.792295Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **<span style = 'color:green'>Interpretable NBeats for Total Sales</span>**<a id ='Total-sale'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\n### **<span style = 'color:brown'>Model Building & Hypertuning of model</span>**<a id ='Building'></a>\n\nLet's take 120 days of horizon and a double of history.","metadata":{}},{"cell_type":"code","source":"back_horizon = 2 * 120\nhorizon = 120\ndate = pd.date_range('2013-01-01', '2017-08-15', freq='D')\ndataset = DataSet(horizon, back_horizon)\ndataset.preprocessing(Total_sales[['sales']].values, Total_sales.index, train_size=0.60, val_size=0.20)","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:25.794296Z","iopub.execute_input":"2022-07-27T08:03:25.794922Z","iopub.status.idle":"2022-07-27T08:03:25.969818Z","shell.execute_reply.started":"2022-07-27T08:03:25.794885Z","shell.execute_reply":"2022-07-27T08:03:25.968854Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Code Credits: [Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}},{"cell_type":"code","source":"def model_builder(hp):\n    \n    hp_n_neurons_trend = hp.Int('neurons_trend', min_value = 4, max_value=16, step=2)\n    hp_n_neurons_seas = hp.Int('neurons_seas', min_value = 4, max_value=16, step=2)\n    hp_periods = hp.Int('periods', min_value = 1, max_value = horizon, step = 30)\n    hp_back_periods = hp.Int('back_periods', min_value = 1, max_value = back_horizon, step = 30)\n    hp_share = hp.Boolean('share')\n    \n    model_nbeats = create_interpretable_nbeats(horizon=horizon, \n                                               back_horizon=back_horizon,\n                                               p_degree=1,   \n                                               trend_n_neurons=hp_n_neurons_trend, \n                                               seasonality_n_neurons=hp_n_neurons_seas, \n                                               periods=hp_periods, \n                                               back_periods=hp_back_periods, \n                                               horizon_fourier_order=hp_periods,\n                                               back_horizon_fourier_order=hp_back_periods,\n                                               n_quantiles=3, \n                                               share=hp_share)\n\n    model_nbeats.compile(loss=QuantileLossError([0.1, 0.5, 0.9]), optimizer=tf.keras.optimizers.Adam(learning_rate=0.001),\n                         metrics=[tf.keras.metrics.MeanAbsoluteError()])\n\n    return model_nbeats\n!pip install -q -U keras-tuner\nimport keras_tuner as kt\ntuner = kt.Hyperband(model_builder,\n                     objective = 'val_loss', \n                     max_epochs = 10,\n                     factor = 3)\n\nimport IPython\nclass ClearTrainingOutput(tf.keras.callbacks.Callback):\n    def on_train_end(*args, **kwargs):\n        IPython.display.clear_output(wait = True)\n\ntuner.search(dataset.X_train, dataset.y_train, epochs = 20, \n             validation_data = (dataset.X_val, dataset.y_val), callbacks = [ClearTrainingOutput()])","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:25.971469Z","iopub.execute_input":"2022-07-27T08:03:25.971994Z","iopub.status.idle":"2022-07-27T08:03:38.601147Z","shell.execute_reply.started":"2022-07-27T08:03:25.971962Z","shell.execute_reply":"2022-07-27T08:03:38.599614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Model Training</span>**<a id ='Training'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\n\nCode Credits: [Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}},{"cell_type":"code","source":"best_hps = tuner.get_best_hyperparameters(num_trials = 1)[0]\n\n# Build the model with the optimal hyperparameters and train it on the data\nmodel = tuner.hypermodel.build(best_hps)\nmodel.build((None, back_horizon))\nmodel.summary()\nfrom tensorflow.keras.callbacks import ReduceLROnPlateau, EarlyStopping\nreduce_lr = ReduceLROnPlateau(monitor='val_loss', factor=0.5, patience=3, min_lr=0.001)\nes_cb = EarlyStopping(monitor='val_loss', min_delta=0,  patience=10, verbose=0, mode='auto', restore_best_weights=True)\n\nhistory = model.fit(dataset.X_train, dataset.y_train, batch_size=64, epochs=20, \n                    callbacks = [reduce_lr, es_cb], validation_data = (dataset.X_val, dataset.y_val))","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:38.603467Z","iopub.execute_input":"2022-07-27T08:03:38.603816Z","iopub.status.idle":"2022-07-27T08:03:47.619257Z","shell.execute_reply.started":"2022-07-27T08:03:38.603784Z","shell.execute_reply":"2022-07-27T08:03:47.618053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plotting training and validation loss during training","metadata":{}},{"cell_type":"code","source":"fig = plt.figure(figsize=(20,5))\n\nax = plt.subplot(131)\n\nepochs = [i for i in range(len(history.history['loss']))]\nax.plot(epochs, history.history['loss'], label='loss')\nax.plot(epochs, history.history['val_loss'], label='val_loss')\nplt.legend()\n\nax = plt.subplot(132)\nax.plot(epochs, history.history['mean_absolute_error'], label='mean_absolute_error')\nax.plot(epochs, history.history['val_mean_absolute_error'], label='val_mean_absolute_error')\nplt.legend()\n\nax = plt.subplot(133)\nax.plot(epochs, history.history['lr'], label='lr')\nb = plt.legend()","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:47.621156Z","iopub.execute_input":"2022-07-27T08:03:47.621521Z","iopub.status.idle":"2022-07-27T08:03:48.147210Z","shell.execute_reply.started":"2022-07-27T08:03:47.621480Z","shell.execute_reply":"2022-07-27T08:03:48.143248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Output Decomposition</span>**<a id ='Decomposition'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\n\nCode Credits: [Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}},{"cell_type":"code","source":"import tensorflow_probability as tfp\ntfd = tfp.distributions\nfrom plotly.subplots import make_subplots\nimport plotly.graph_objects as go\noutputs = []\nseasonality = []\ntrend = []\nfor sample in range(30):\n    results = model(dataset.X_test)\n    outputs.append(results)\n    seasonality.append(model.seasonality)\n    trend.append(model.trend)\n    \nplot_results_nbeats(y_pred=tf.stack(outputs),\n                    date_outputs=dataset.test_date,\n                    y_true=dataset.y_test,\n                    date_history=None,\n                    seasonality=tf.stack(seasonality),\n                    trend=tf.stack(trend))","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:48.149526Z","iopub.execute_input":"2022-07-27T08:03:48.150007Z","iopub.status.idle":"2022-07-27T08:03:49.105885Z","shell.execute_reply.started":"2022-07-27T08:03:48.149959Z","shell.execute_reply":"2022-07-27T08:03:49.104771Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Model Forecast</span>**<a id ='Forecast'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)","metadata":{}},{"cell_type":"code","source":"outputs = []\nseasonality = []\ntrend = []\nfor sample in range(50):\n    results = model(Total_sales[['sales']].iloc[-240:].T.values)\n    outputs.append(results)\n    seasonality.append(model.seasonality)\n    trend.append(model.trend)","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:49.107289Z","iopub.execute_input":"2022-07-27T08:03:49.107625Z","iopub.status.idle":"2022-07-27T08:03:50.035371Z","shell.execute_reply.started":"2022-07-27T08:03:49.107592Z","shell.execute_reply":"2022-07-27T08:03:50.034087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"75 days forecast from the end date of train data. ","metadata":{}},{"cell_type":"code","source":"date_forecasting = np.expand_dims(pd.date_range('2017-08-16', '2017-10-31', freq= 'D'), axis=0)\n\nplot_results_nbeats(y_pred=tf.stack(outputs),\n                    date_outputs=date_forecasting,\n                    y_true=Total_sales[['sales']].T.values,\n                    #y_true=df.T.values,\n                    date_history=date,\n                    seasonality=tf.stack(seasonality),\n                    trend=tf.stack(trend))","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:50.037700Z","iopub.execute_input":"2022-07-27T08:03:50.038243Z","iopub.status.idle":"2022-07-27T08:03:50.215444Z","shell.execute_reply.started":"2022-07-27T08:03:50.038189Z","shell.execute_reply":"2022-07-27T08:03:50.214213Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **<span style = 'color:green'>Interpretable NBeats for Daily Store Sales Data</span>**<a id ='Store-sale'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)\n\n### **<span style = 'color:brown'>Model Building & Training</span>**<a id ='Building2'></a>\n\nRepeat the steps above to build, hypertune, train & use the model for prediction.","metadata":{}},{"cell_type":"code","source":"back_horizon = 2 * 120\nhorizon = 120\ndate = pd.date_range('2013-01-01', '2017-08-15', freq='D')\ndataset = DataSet(horizon, back_horizon)\ndataset.preprocessing(store_sales.values, store_sales.index, train_size=0.60, val_size=0.20)","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:03:50.217415Z","iopub.execute_input":"2022-07-27T08:03:50.218174Z","iopub.status.idle":"2022-07-27T08:04:00.411434Z","shell.execute_reply.started":"2022-07-27T08:03:50.218123Z","shell.execute_reply":"2022-07-27T08:04:00.410090Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def model_builder(hp):\n    \n    hp_n_neurons_trend = hp.Int('neurons_trend', min_value = 4, max_value=16, step=2)\n    hp_n_neurons_seas = hp.Int('neurons_seas', min_value = 4, max_value=16, step=2)\n    hp_periods = hp.Int('periods', min_value = 1, max_value = horizon, step = 30)\n    hp_back_periods = hp.Int('back_periods', min_value = 1, max_value = back_horizon, step = 30)\n    hp_share = hp.Boolean('share')\n    \n    model_nbeats = create_interpretable_nbeats(horizon=horizon, \n                                               back_horizon=back_horizon,\n                                               p_degree=1,   \n                                               trend_n_neurons=hp_n_neurons_trend, \n                                               seasonality_n_neurons=hp_n_neurons_seas, \n                                               periods=hp_periods, \n                                               back_periods=hp_back_periods, \n                                               horizon_fourier_order=hp_periods,\n                                               back_horizon_fourier_order=hp_back_periods,\n                                               n_quantiles=3, \n                                               share=hp_share)\n\n    model_nbeats.compile(loss=QuantileLossError([0.1, 0.5, 0.9]), optimizer=tf.keras.optimizers.Adam(learning_rate=0.001),\n                         metrics=[tf.keras.metrics.MeanAbsoluteError()])\n\n    return model_nbeats\n!pip install -q -U keras-tuner\nimport keras_tuner as kt\ntuner = kt.Hyperband(model_builder,\n                     objective = 'val_loss', \n                     max_epochs = 10,\n                     factor = 3)\n\nimport IPython\nclass ClearTrainingOutput(tf.keras.callbacks.Callback):\n    def on_train_end(*args, **kwargs):\n        IPython.display.clear_output(wait = True)\n\ntuner.search(dataset.X_train, dataset.y_train, epochs = 20, \n             validation_data = (dataset.X_val, dataset.y_val), callbacks = [ClearTrainingOutput()])","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:04:00.415472Z","iopub.execute_input":"2022-07-27T08:04:00.415851Z","iopub.status.idle":"2022-07-27T08:04:12.834524Z","shell.execute_reply.started":"2022-07-27T08:04:00.415817Z","shell.execute_reply":"2022-07-27T08:04:12.833018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"best_hps = tuner.get_best_hyperparameters(num_trials = 1)[0]\n\n# Build the model with the optimal hyperparameters and train it on the data\nmodel = tuner.hypermodel.build(best_hps)\nmodel.build((None, back_horizon))\nmodel.summary()\nfrom tensorflow.keras.callbacks import ReduceLROnPlateau, EarlyStopping\nreduce_lr = ReduceLROnPlateau(monitor='val_loss', factor=0.5, patience=3, min_lr=0.001)\nes_cb = EarlyStopping(monitor='val_loss', min_delta=0,  patience=10, verbose=0, mode='auto', restore_best_weights=True)\n\nhistory = model.fit(dataset.X_train, dataset.y_train, batch_size=64, epochs=20, \n                    callbacks = [reduce_lr, es_cb], validation_data = (dataset.X_val, dataset.y_val))","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:04:12.836648Z","iopub.execute_input":"2022-07-27T08:04:12.837034Z","iopub.status.idle":"2022-07-27T08:06:10.622868Z","shell.execute_reply.started":"2022-07-27T08:04:12.836998Z","shell.execute_reply":"2022-07-27T08:06:10.621496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(20,5))\n\nax = plt.subplot(131)\n\nepochs = [i for i in range(len(history.history['loss']))]\nax.plot(epochs, history.history['loss'], label='loss')\nax.plot(epochs, history.history['val_loss'], label='val_loss')\nplt.legend()\n\nax = plt.subplot(132)\nax.plot(epochs, history.history['mean_absolute_error'], label='mean_absolute_error')\nax.plot(epochs, history.history['val_mean_absolute_error'], label='val_mean_absolute_error')\nplt.legend()\n\nax = plt.subplot(133)\nax.plot(epochs, history.history['lr'], label='lr')\nb = plt.legend()","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:06:10.626923Z","iopub.execute_input":"2022-07-27T08:06:10.628456Z","iopub.status.idle":"2022-07-27T08:06:11.101193Z","shell.execute_reply.started":"2022-07-27T08:06:10.628410Z","shell.execute_reply":"2022-07-27T08:06:11.099397Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import tensorflow_probability as tfp\ntfd = tfp.distributions\nfrom plotly.subplots import make_subplots\nimport plotly.graph_objects as go\noutputs = []\nseasonality = []\ntrend = []\nfor sample in range(30):\n    results = model(dataset.X_test)\n    outputs.append(results)\n    seasonality.append(model.seasonality)\n    trend.append(model.trend)\n    \nplot_results_nbeats(y_pred=tf.stack(outputs),\n                    date_outputs=dataset.test_date,\n                    y_true=dataset.y_test,\n                    date_history=None,\n                    seasonality=tf.stack(seasonality),\n                    trend=tf.stack(trend))","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:06:11.103213Z","iopub.execute_input":"2022-07-27T08:06:11.103910Z","iopub.status.idle":"2022-07-27T08:06:19.199764Z","shell.execute_reply.started":"2022-07-27T08:06:11.103857Z","shell.execute_reply":"2022-07-27T08:06:19.198582Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### **<span style = 'color:brown'>Model Forecast of Daily Store Sales</span>**<a id ='Forecast2'></a>\n[<div style=\"text-align: right\"> Back to Table of contents</div>](#Table)","metadata":{}},{"cell_type":"code","source":"outputs = []\nseasonality = []\ntrend = []\nfor sample in range(50):\n    results = model(store_sales.iloc[-240:].T.values)\n    outputs.append(results)\n    seasonality.append(model.seasonality)\n    trend.append(model.trend)","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:06:19.201685Z","iopub.execute_input":"2022-07-27T08:06:19.202415Z","iopub.status.idle":"2022-07-27T08:06:20.171032Z","shell.execute_reply.started":"2022-07-27T08:06:19.202367Z","shell.execute_reply":"2022-07-27T08:06:20.169815Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"date_forecasting = np.expand_dims(pd.date_range('2017-08-16', '2017-10-31', freq= 'D'), axis=0)\n\nplot_results_nbeats(y_pred=tf.stack(outputs),\n                    date_outputs=date_forecasting,\n                    y_true=store_sales.T.values,\n                    #y_true=df.T.values,\n                    date_history=date,\n                    seasonality=tf.stack(seasonality),\n                    trend=tf.stack(trend))","metadata":{"execution":{"iopub.status.busy":"2022-07-27T08:06:20.172707Z","iopub.execute_input":"2022-07-27T08:06:20.173019Z","iopub.status.idle":"2022-07-27T08:06:20.356112Z","shell.execute_reply.started":"2022-07-27T08:06:20.172987Z","shell.execute_reply":"2022-07-27T08:06:20.355053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **Reference:**\n[Forecast with N-BEATS || Interpretable model](https://www.kaggle.com/code/gatandubuc/forecast-with-n-beats-interpretable-model/notebook)","metadata":{}}]}