{"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":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-07-23T09:06:33.346790Z","iopub.execute_input":"2022-07-23T09:06:33.347415Z","iopub.status.idle":"2022-07-23T09:06:33.384458Z","shell.execute_reply.started":"2022-07-23T09:06:33.347323Z","shell.execute_reply":"2022-07-23T09:06:33.383522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predict Future Sales\n\n1. Import Data and Data Preprocessing\n\n2. Exploratory Data Analysis\n\n3. Data Model and Prediction","metadata":{}},{"cell_type":"code","source":"# Import Sales Train as Main Table\nMain_df = pd.read_csv('/kaggle/input/competitive-data-science-predict-future-sales/sales_train.csv')\nMain_df['date'] = pd.to_datetime(Main_df['date'], format = '%d.%m.%Y')\n\n# Create Amount\nMain_df['amount'] = Main_df['item_price'] * Main_df['item_cnt_day']\nMain_df.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:06:33.386083Z","iopub.execute_input":"2022-07-23T09:06:33.386579Z","iopub.status.idle":"2022-07-23T09:06:36.658484Z","shell.execute_reply.started":"2022-07-23T09:06:33.386533Z","shell.execute_reply":"2022-07-23T09:06:36.657672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Data Checking\nprint('Shape of Sales Data:', Main_df.shape)\nprint('Missing Entries:',Main_df.isna().sum().sum())\nprint('Data Type of Sales Data:\\n', Main_df.dtypes, '\\n')","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:06:36.659609Z","iopub.execute_input":"2022-07-23T09:06:36.659833Z","iopub.status.idle":"2022-07-23T09:06:36.703900Z","shell.execute_reply.started":"2022-07-23T09:06:36.659808Z","shell.execute_reply":"2022-07-23T09:06:36.702086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Import Supplment Information\nShop_df = pd.read_csv('/kaggle/input/competitive-data-science-predict-future-sales/shops.csv')\nShop_df.shape","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:06:36.705740Z","iopub.execute_input":"2022-07-23T09:06:36.706051Z","iopub.status.idle":"2022-07-23T09:06:36.718495Z","shell.execute_reply.started":"2022-07-23T09:06:36.706019Z","shell.execute_reply":"2022-07-23T09:06:36.717624Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Category_df = pd.read_csv('/kaggle/input/competitive-data-science-predict-future-sales/item_categories.csv')\nCategory_df.shape","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:06:36.719957Z","iopub.execute_input":"2022-07-23T09:06:36.720279Z","iopub.status.idle":"2022-07-23T09:06:36.732102Z","shell.execute_reply.started":"2022-07-23T09:06:36.720237Z","shell.execute_reply":"2022-07-23T09:06:36.731348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Item_df = pd.read_csv('/kaggle/input/competitive-data-science-predict-future-sales/items.csv')\nItem_df.shape","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:06:36.733729Z","iopub.execute_input":"2022-07-23T09:06:36.733963Z","iopub.status.idle":"2022-07-23T09:06:36.809854Z","shell.execute_reply.started":"2022-07-23T09:06:36.733935Z","shell.execute_reply":"2022-07-23T09:06:36.808947Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Merge the Supplemntory Information into the Main Table\nprint('Before Merging:', Main_df.shape)\nMain_df = pd.merge(Main_df, Shop_df, on = 'shop_id', how = 'left')\nMain_df = pd.merge(Main_df, Item_df, on = 'item_id', how = 'left')\nMain_df = pd.merge(Main_df, Category_df, on = 'item_category_id', how = 'left')\nprint('After Merging:', Main_df.shape)","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:06:36.810950Z","iopub.execute_input":"2022-07-23T09:06:36.811548Z","iopub.status.idle":"2022-07-23T09:06:38.327745Z","shell.execute_reply.started":"2022-07-23T09:06:36.811513Z","shell.execute_reply":"2022-07-23T09:06:38.326719Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Sum the Negative Count First\nMain_df['Year Month'] = Main_df['date'].apply(lambda x:x.strftime('%Y%m'))\nMain_df = Main_df.groupby(['Year Month',\n                           'shop_id',\n                           'item_category_id',\n                           'item_id'],as_index = False).agg({'item_cnt_day':'sum',\n                                                             'amount':'sum'})","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:06:38.329265Z","iopub.execute_input":"2022-07-23T09:06:38.329491Z","iopub.status.idle":"2022-07-23T09:07:08.706224Z","shell.execute_reply.started":"2022-07-23T09:06:38.329464Z","shell.execute_reply":"2022-07-23T09:07:08.705153Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Exploratory Data Analysis","metadata":{}},{"cell_type":"code","source":"# TSA from Statsmodels\nimport statsmodels.api as sm\nimport statsmodels.formula.api as smf\nimport statsmodels.tsa.api as smt\n\n# Display and Plotting\nimport matplotlib.pylab as plt\nimport seaborn as sns\n\nsns.set_theme()","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:07:08.707890Z","iopub.execute_input":"2022-07-23T09:07:08.708254Z","iopub.status.idle":"2022-07-23T09:07:10.576190Z","shell.execute_reply.started":"2022-07-23T09:07:08.708211Z","shell.execute_reply":"2022-07-23T09:07:10.575165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15, 7))\nsns.lineplot(data = Main_df, x = 'Year Month', y =  'item_cnt_day').set_title('Items Sold by Date')\nplt.xticks(rotation = 90)\nplt.show()\n\nplt.figure(figsize=(15, 7))\nsns.lineplot(data = Main_df, x = 'Year Month', y = 'amount').set_title('Amount of Sales by Date')\nplt.xticks(rotation = 90)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:07:10.579072Z","iopub.execute_input":"2022-07-23T09:07:10.580063Z","iopub.status.idle":"2022-07-23T09:07:50.995818Z","shell.execute_reply.started":"2022-07-23T09:07:10.580020Z","shell.execute_reply":"2022-07-23T09:07:50.994753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Findings:\n\n1. The graph of quantity and amount demonstrated a similar seasonal pattern, where sales are mainly flat across the year, except an peak in December.\n2. In general, the quantity sold has demonstrated a minor decrease, especially in 2015. However, the sales amount shows a rise in opposite.\n3. In 2015, the quantity seems more volatile comparing to previous 2 years.","metadata":{}},{"cell_type":"markdown","source":"# Data Preprocessing","metadata":{}},{"cell_type":"markdown","source":"# Dimensonal Reduction","metadata":{}},{"cell_type":"markdown","source":"As this dataset is time-series data, the internal sequences and permutations matter. The vector cannot be broken into separate rows like normal regression analysis, where the data will lose its serial feature. Given this circumstance, an autoregression model will suit our purpose well.\n\nHowever, the prediction requires not only the prediction of item_cnt_date in Nov 2015, but it also considers the differentiation between both shops and products. This requirement greatly increased the dimensionality complexity. Under the massive combinations, we should consider the resource allocation and time complexity beneath.\n\nTo accomplish the goal, we will employ a generative model: first conduct a dimensional reduction to differentiate some important clusters, upon which multiple autoregression models will be built for our prediction.","metadata":{}},{"cell_type":"markdown","source":"As the prediction interval is by month, we do not have to consider day differentiation. Let us aggregate the data into year month level for our cluster analysis.","metadata":{}},{"cell_type":"code","source":"# Pivot the Data\nYear_Month = Main_df['Year Month'].unique()\nMain_df = Main_df.groupby(['Year Month','shop_id','item_category_id','item_id']).agg({'item_cnt_day': 'nunique'}).unstack('Year Month').reset_index()\nMain_df.columns = np.concatenate((['shop_id','item_category_id','item_id'],Year_Month.tolist()))\nMain_df","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:07:50.997075Z","iopub.execute_input":"2022-07-23T09:07:50.997326Z","iopub.status.idle":"2022-07-23T09:07:53.170269Z","shell.execute_reply.started":"2022-07-23T09:07:50.997298Z","shell.execute_reply":"2022-07-23T09:07:53.169376Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Missing Data","metadata":{}},{"cell_type":"code","source":"# How Many Column Missing?\nprint(Main_df.isna().sum()/ len(Main_df))","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:07:53.171502Z","iopub.execute_input":"2022-07-23T09:07:53.171889Z","iopub.status.idle":"2022-07-23T09:07:53.211264Z","shell.execute_reply.started":"2022-07-23T09:07:53.171853Z","shell.execute_reply":"2022-07-23T09:07:53.210389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Missing values are common in our dataset. To handle the sparse matrix, we must consider our goals and data features. Here are some questions when coming up with a solution.\n\n1. Will the removal of data help, or pose great impacts on our analysis?\n2. If we cannot remove the data, what will be the best way of substitution for our data feature?\n\nMy thought is the removal of the missing data is not a good choice. We have two ways of removals:\n\nRemoving by row or column will either undermine 1) the prediction by shops and products and 2) break the time series dependence.\n\nThen what's next is the substitution. There are three major solutions, replacing them with 0, average, or previous values.\n\n1. replacing the average may not be a good option, as it is referring to the mean of a time series which may undermine the sequential feature.\n2. replacing with previous values sounds more rational as it considers the sequential effect, but we are not sure about the daily patterns involved and we may create many additional sales that undermine our model performance.\n3. replacing with 0 seems a good choice, as it just keeps the original patterns, and it is common to have no sale on a single day, especially under massive combinations\n\nWe admit that there are trade-offs for each method. And we do not have the relevant domain knowledge so we can just do it in a more conservative way here.","metadata":{}},{"cell_type":"code","source":"# Replace with 0\nMain_df = Main_df.replace(np.nan, 0)","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:07:53.212389Z","iopub.execute_input":"2022-07-23T09:07:53.212608Z","iopub.status.idle":"2022-07-23T09:07:53.297679Z","shell.execute_reply.started":"2022-07-23T09:07:53.212582Z","shell.execute_reply":"2022-07-23T09:07:53.296749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let use an elbow method to evaluate the optimial k cluster for the prediction.","metadata":{}},{"cell_type":"code","source":"from sklearn.cluster import KMeans\n\nk = [i for i in range(1, 11)]\ndistance = []\n\nfor i in k:\n    kmeans = KMeans(n_clusters=i, random_state=0).fit(Main_df.iloc[:,3:])\n    distance.append(kmeans.inertia_)\n    \nplt.title('Elbow Method of K-Means Clustering')\nplt.plot(k,\n         distance\n        )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:07:53.299177Z","iopub.execute_input":"2022-07-23T09:07:53.299571Z","iopub.status.idle":"2022-07-23T09:10:25.563568Z","shell.execute_reply.started":"2022-07-23T09:07:53.299526Z","shell.execute_reply":"2022-07-23T09:10:25.562796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The elbow method shows that the critical point for optimal k is around 4, where additional clusters do not help to reduce the variances. We will then map the cluster for each combination, and create models upon these 3 clusters generated by the k-mean cluster model.","metadata":{}},{"cell_type":"code","source":"kmeans = KMeans(n_clusters=4, random_state=0).fit(Main_df.iloc[:,3:])\nMain_df['Cluster'] = kmeans.labels_","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:10:25.565321Z","iopub.execute_input":"2022-07-23T09:10:25.565817Z","iopub.status.idle":"2022-07-23T09:10:37.125452Z","shell.execute_reply.started":"2022-07-23T09:10:25.565780Z","shell.execute_reply":"2022-07-23T09:10:37.124537Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Plot_df = Main_df.iloc[:,3:].groupby('Cluster',as_index = False).sum().melt('Cluster')\nPlot_df.columns = ['Cluster','Year Month','item_cnt_day']\n\nplt.figure(figsize=(15, 7))\nsns.lineplot(data = Plot_df, x = 'Year Month', y = 'item_cnt_day', hue = 'Cluster').set_title('Total item_cnt_day by Cluster')\nplt.xticks(rotation = 90)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:10:37.126589Z","iopub.execute_input":"2022-07-23T09:10:37.126936Z","iopub.status.idle":"2022-07-23T09:10:37.826263Z","shell.execute_reply.started":"2022-07-23T09:10:37.126902Z","shell.execute_reply":"2022-07-23T09:10:37.825570Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Show Distribution\nMain_df['Cluster'].value_counts().sort_index()","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:10:37.827431Z","iopub.execute_input":"2022-07-23T09:10:37.827699Z","iopub.status.idle":"2022-07-23T09:10:37.838947Z","shell.execute_reply.started":"2022-07-23T09:10:37.827655Z","shell.execute_reply":"2022-07-23T09:10:37.837894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The cluster plots show some differentiation between clusters:\n\n1. Two clusters share similar patterns where a peak appeared in Dec 2013, while one cluster had a strong re-bounce since June 2015.\n2. A cluster shows a gradual downtrend since 2013.\n3. A cluster shows a strong surge in 2014.\n\nThen the next step is to build the corresponding model upon the above clusters.\n(We did not include specific cluster name and index as it may vary due to random assigning)","metadata":{}},{"cell_type":"markdown","source":"# Time Series Model","metadata":{}},{"cell_type":"markdown","source":"Let us build forecastting upon the above clustering results!","metadata":{}},{"cell_type":"code","source":"# Preprocessing \n\nModel_df = Main_df.iloc[:,3:].melt('Cluster').groupby(['Cluster','variable']).mean().unstack('Cluster').reset_index()\nModel_df.columns = ['Year Month','Cluster 0','Cluster 1','Cluster 2','Cluster 3']\nModel_df['Year Month'] = pd.date_range('20130101','20151031', freq = 'm')\nModel_df = Model_df.set_index('Year Month')","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:10:37.840002Z","iopub.execute_input":"2022-07-23T09:10:37.840204Z","iopub.status.idle":"2022-07-23T09:10:40.420547Z","shell.execute_reply.started":"2022-07-23T09:10:37.840178Z","shell.execute_reply":"2022-07-23T09:10:40.419551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot Summary Graph:\n\ndef tsplot(y, lags=None, title='', figsize=(14, 8)):\n    '''Examine the patterns of ACF and PACF, along with the time series plot and histogram.\n    \n    Original source: https://tomaugspurger.github.io/modern-7-timeseries.html\n    '''\n    fig = plt.figure(figsize=figsize)\n    layout = (2, 2)\n    ts_ax   = plt.subplot2grid(layout, (0, 0))\n    hist_ax = plt.subplot2grid(layout, (0, 1))\n    acf_ax  = plt.subplot2grid(layout, (1, 0))\n    pacf_ax = plt.subplot2grid(layout, (1, 1))\n    \n    y.plot(ax=ts_ax)\n    ts_ax.set_title(title)\n    y.plot(ax=hist_ax, kind='hist', bins=25)\n    hist_ax.set_title('Histogram')\n    smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)\n    smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)\n    [ax.set_xlim(0) for ax in [acf_ax, pacf_ax]]\n    sns.despine()\n    plt.tight_layout()\n    return ts_ax, acf_ax, pacf_ax","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:10:40.422158Z","iopub.execute_input":"2022-07-23T09:10:40.422709Z","iopub.status.idle":"2022-07-23T09:10:40.434011Z","shell.execute_reply.started":"2022-07-23T09:10:40.422644Z","shell.execute_reply":"2022-07-23T09:10:40.433087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in Model_df.columns:\n    tsplot(Model_df[i], lags=12, title=i, figsize=(14, 8))","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:10:40.435618Z","iopub.execute_input":"2022-07-23T09:10:40.436534Z","iopub.status.idle":"2022-07-23T09:10:46.152907Z","shell.execute_reply.started":"2022-07-23T09:10:40.436482Z","shell.execute_reply":"2022-07-23T09:10:46.152008Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Train Test Split (9:1)\nTrain_x = Model_df[:round(len(Model_df)*0.9)]\nTest_x = Model_df[round(len(Model_df)*0.1):]","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:10:46.154473Z","iopub.execute_input":"2022-07-23T09:10:46.154823Z","iopub.status.idle":"2022-07-23T09:10:46.163419Z","shell.execute_reply.started":"2022-07-23T09:10:46.154779Z","shell.execute_reply":"2022-07-23T09:10:46.162339Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Model Estimation\n\n# Fit the model\nfor i in Train_x.columns:\n    arima = sm.tsa.arima.ARIMA(Train_x[i], order=(4,0,2),  freq = 'M')\n    model_results = arima.fit()\n    print(model_results.summary())","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:10:46.164584Z","iopub.execute_input":"2022-07-23T09:10:46.164850Z","iopub.status.idle":"2022-07-23T09:10:48.194293Z","shell.execute_reply.started":"2022-07-23T09:10:46.164815Z","shell.execute_reply":"2022-07-23T09:10:48.193345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Opps, we observe that most of the default parameters are not working well. Most of the autoregression and moving average do not show a significant p-value (p-value < 0.05). We cannot accept the hypotheses that the autoregression is highly significant to the next-step prediction. Let us try to tune the hyperparameters by AIC and BIC as our evaluation metrics!","metadata":{}},{"cell_type":"code","source":"# Alternative model selection method, limited to only searching AR and MA parameters\n\naic_order = []\nbic_order = []\n\nfor i in Model_df.columns:\n    \n    train_results = sm.tsa.arma_order_select_ic(Train_x[i], ic=['aic', 'bic'], trend='n', max_ar=5, max_ma=5)\n    aic_order.append(train_results.aic_min_order)\n    bic_order.append(train_results.bic_min_order)\n\n    print(f'{i} - AIC', train_results.aic_min_order)\n    print(f'{i} - BIC', train_results.bic_min_order)","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:10:48.196427Z","iopub.execute_input":"2022-07-23T09:10:48.197201Z","iopub.status.idle":"2022-07-23T09:11:33.806731Z","shell.execute_reply.started":"2022-07-23T09:10:48.197146Z","shell.execute_reply":"2022-07-23T09:11:33.805732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Summarized the AIC Order for Each Model\npd.DataFrame({'Cluster': Model_df.columns,\n              'AIC Order': aic_order,\n              'BIC Order': bic_order})","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:11:33.812677Z","iopub.execute_input":"2022-07-23T09:11:33.816786Z","iopub.status.idle":"2022-07-23T09:11:33.845815Z","shell.execute_reply.started":"2022-07-23T09:11:33.816707Z","shell.execute_reply":"2022-07-23T09:11:33.844706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Model Estimation\n\n# Fit the model\nfor i in range(4):\n    arima = sm.tsa.arima.ARIMA(Train_x.iloc[:, i], order=(bic_order[i][0],\n                                                          0,\n                                                          bic_order[i][1]),  freq = 'M')\n    model_results = arima.fit()\n    print(model_results.summary())\n    \n    # Re-run the above statistical tests, and more. To be used when selecting viable models.\n\n    het_method='breakvar'\n    norm_method='jarquebera'\n    sercor_method='ljungbox'\n\n    (het_stat, het_p) = model_results.test_heteroskedasticity(het_method)[0]\n    norm_stat, norm_p, skew, kurtosis = model_results.test_normality(norm_method)[0]\n    sercor_stat, sercor_p = model_results.test_serial_correlation(method=sercor_method)[0]\n    sercor_stat = sercor_stat[-1] # last number for the largest lag\n    sercor_p = sercor_p[-1] # last number for the largest lag\n\n    # Run Durbin-Watson test on the standardized residuals.\n    # The statistic is approximately equal to 2*(1-r), where r is the sample autocorrelation of the residuals.\n    # Thus, for r == 0, indicating no serial correlation, the test statistic equals 2.\n    # This statistic will always be between 0 and 4. The closer to 0 the statistic,\n    # the more evidence for positive serial correlation. The closer to 4,\n    # the more evidence for negative serial correlation.\n    # Essentially, below 1 or above 3 is bad.\n    dw = sm.stats.stattools.durbin_watson(model_results.filter_results.standardized_forecasts_error[0, model_results.loglikelihood_burn:])\n\n    # check whether roots are outside the unit circle (we want them to be);\n    # will be True when AR is not used (i.e., AR order = 0)\n    arroots_outside_unit_circle = np.all(np.abs(model_results.arroots) > 1)\n    # will be True when MA is not used (i.e., MA order = 0)\n    maroots_outside_unit_circle = np.all(np.abs(model_results.maroots) > 1)\n\n    print('Test heteroskedasticity of residuals ({}): stat={:.3f}, p={:.3f}'.format(het_method, het_stat, het_p));\n    print('\\nTest normality of residuals ({}): stat={:.3f}, p={:.3f}'.format(norm_method, norm_stat, norm_p));\n    print('\\nTest serial correlation of residuals ({}): stat={:.3f}, p={:.3f}'.format(sercor_method, sercor_stat, sercor_p));\n    print('\\nDurbin-Watson test on residuals: d={:.2f}\\n\\t(NB: 2 means no serial correlation, 0=pos, 4=neg)'.format(dw))\n    print('\\nTest for all AR roots outside unit circle (>1): {}'.format(arroots_outside_unit_circle))\n    print('\\nTest for all MA roots outside unit circle (>1): {}'.format(maroots_outside_unit_circle))\n\n    model_results.plot_diagnostics(figsize=(16, 12))\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:11:33.853102Z","iopub.execute_input":"2022-07-23T09:11:33.856508Z","iopub.status.idle":"2022-07-23T09:11:38.781603Z","shell.execute_reply.started":"2022-07-23T09:11:33.856430Z","shell.execute_reply":"2022-07-23T09:11:38.780969Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_rmse(y, y_hat):\n    '''Root Mean Square Error\n    '''\n    mse = np.mean((y - y_hat)**2)\n    return np.sqrt(mse)\n\ndef get_mape(y, y_hat):\n    '''Mean Absolute Percent Error\n    '''\n    perc_err = (100*(y - y_hat))/y\n    return np.mean(abs(perc_err))\n\ndef get_mase(y, y_hat):\n    '''Mean Absolute Scaled Error\n    '''\n    abs_err = abs(y - y_hat)\n    dsum=sum(abs(y[1:] - y_hat[1:]))\n    t = len(y)\n    denom = (1/(t - 1))* dsum\n    return np.mean(abs_err/denom)","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:11:38.782605Z","iopub.execute_input":"2022-07-23T09:11:38.783292Z","iopub.status.idle":"2022-07-23T09:11:38.790874Z","shell.execute_reply.started":"2022-07-23T09:11:38.783257Z","shell.execute_reply":"2022-07-23T09:11:38.789861Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prediction = []\n\nfor i in range(4):\n    arima = sm.tsa.arima.ARIMA(Train_x.iloc[:, i], order=(bic_order[i][0],\n                                                          0,\n                                                          bic_order[i][1]),  freq = 'M')\n    model_results = arima.fit()\n    \n    fig, ax1 = plt.subplots(nrows=1, ncols=1, figsize=(12, 8))\n    \n    ax1.plot(Train_x.iloc[:, i], label='In-sample data', linestyle='-')\n    # subtract 1 only to connect it to previous point in the graph\n    ax1.plot(Test_x.iloc[:, i], label='Held-out data', linestyle='--')\n\n    # yes DatetimeIndex\n    pred_begin = Train_x.iloc[:, i].index[model_results.loglikelihood_burn]\n    pred_end = Test_x.iloc[:, i].index[-1]\n    pred = model_results.get_prediction(start = 0, end = len(Model_df))\n\n    pred_mean = pred.predicted_mean[1:]\n    pred_mean.index = Model_df.index\n    pred_ci = pred.conf_int(alpha=0.05)[1:]\n    pred_ci.index = Model_df.index\n    \n    prediction.append(pred_mean.to_list()[-1])\n    \n    print('RMSE:', get_rmse(Model_df.iloc[:, i], pred_mean))\n    print('MAPE:', get_mape(Model_df.iloc[:, i], pred_mean))\n    print('MASE:', get_mase(Model_df.iloc[:, i], pred_mean))\n    \n    ax1.plot(pred_mean[1:], 'r', alpha=.6, label='Predicted values')\n    ax1.fill_between(pred_ci[1:].index,\n                     pred_ci.iloc[1:, 0],\n                     pred_ci.iloc[1:, 1], color='k', alpha=.2)\n\n    ax1.legend(loc='best');\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:11:38.791898Z","iopub.execute_input":"2022-07-23T09:11:38.792117Z","iopub.status.idle":"2022-07-23T09:11:40.816389Z","shell.execute_reply.started":"2022-07-23T09:11:38.792091Z","shell.execute_reply":"2022-07-23T09:11:40.815542Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"After tuning the hyperparameters, the p-value is now significant and the autoregression as 1 is quite relevant to the value in the next time step.\n\n\nThe first two clusters are working quite good on predicting the next time step, but the last 2 are not working fine. It may be caused by the homogeneity problem within the cluster, where the first two are quite homogeneous while the last 2 are quite heterogeneous.","metadata":{}},{"cell_type":"code","source":"# Predict for Submission\nPrediction_by_Cluster = pd.DataFrame({'Cluster': [i for i in range(4)],\n                                      'Prediction for Nov 2015':prediction})\nPrediction_by_Cluster.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:11:40.819606Z","iopub.execute_input":"2022-07-23T09:11:40.819860Z","iopub.status.idle":"2022-07-23T09:11:40.830733Z","shell.execute_reply.started":"2022-07-23T09:11:40.819831Z","shell.execute_reply":"2022-07-23T09:11:40.829874Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As you may be aware, the prediction seems abnormal as they include decimal places. This is a drawback when using a generative model to make an aggregated prediction. It will flatten and reduce all predictions to the same under the same cluster. It may create illogical predictions like here, where product counts should be integers. We also found that the clustering is difficult to tackle cold start problem and sparse matrix. Our approach turns out to seem not very feasible, but it is still a good try to consider dimensionality problems and space-time complexity during time-series modelling!\n\nIt is also the reality of retail industry has thousands of products, which makes the prediction and computation difficult. There is some Python package that aims to make the prediction more efficient and scalable, e.g. Prophet. Apart from the package, another issue matters are the use case, whether it is for general business planning or reporting, or it is a real-world inventory management problem. We can scale down like the generative model here, but the latter involves down to product level but not very much about statistical inferences and machine learning explainability, where ML flow and MLops may be helpful to regularize the prediction for the purpose! These are some considerations for our to work on in the future and under a real-world environment!\n","metadata":{}},{"cell_type":"code","source":"# Merge \nSubmission = pd.read_csv('../input/competitive-data-science-predict-future-sales/test.csv')\nSubmission = pd.merge(Submission, Main_df[['shop_id','item_id','Cluster']], on = ['shop_id','item_id'], how = 'left')\n\n# Map with Cluster\nSubmission = pd.merge(Submission, Prediction_by_Cluster, on = 'Cluster', how = 'left')\n\n# If Cold Start, then fill 0\nSubmission['Prediction for Nov 2015'] = Submission['Prediction for Nov 2015'].fillna(0)\nSubmission['item_cnt_month'] = Submission['Prediction for Nov 2015']\n\nSubmission = Submission[['ID','item_cnt_month']]\nSubmission","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:21:00.484061Z","iopub.execute_input":"2022-07-23T09:21:00.484455Z","iopub.status.idle":"2022-07-23T09:21:00.712507Z","shell.execute_reply.started":"2022-07-23T09:21:00.484418Z","shell.execute_reply":"2022-07-23T09:21:00.711903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Submission.to_csv('submission.csv',index= False)","metadata":{"execution":{"iopub.status.busy":"2022-07-23T09:21:09.376957Z","iopub.execute_input":"2022-07-23T09:21:09.377511Z","iopub.status.idle":"2022-07-23T09:21:09.979903Z","shell.execute_reply.started":"2022-07-23T09:21:09.377464Z","shell.execute_reply":"2022-07-23T09:21:09.978949Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Reference:\n\nJeffery Yau at https://www.youtube.com/watch?v=_vQ0W_qXMxk. ","metadata":{}}]}