<jupyter_start><jupyter_text>World Wide Products Inc.Shipping and delivering to a place near you1. Set up a data science project structure in a new git repository in your GitHub account2. Download the product demand data set from https://www.kaggle.com/felixzhao/productdemandforecasting3. Load the data set into panda data frames4. Formulate one or two ideas on how feature engineering would help the data set to establish additional value using exploratory data analysis5. Build one or more forecasting models to determine the demand for a particular product using the other columns as features6. Document your process and results7. Commit your notebook, source code, visualizations and other supporting files to the git repository in GitHub<jupyter_code>%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import datetime

from scipy.stats import norm
import statsmodels.api as sm
import matplotlib.pyplot as plt

import requests
from io import BytesIO
from zipfile import ZipFile
from sklearn import preprocessing

import warnings
from sklearn.metrics import mean_squared_error

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
from sklearn.preprocessing import MinMaxScaler

from fbprophet import Prophet
productDemandDf = pd.read_csv('../data/productdemandforecasting.zip')
productDemandDf.tail(5)
productDemandDf['Product_Code'].nunique()# seeing the number of unique products
productDemandDf.dtypes<jupyter_output><empty_output><jupyter_text>**Since the date feature is in generic object datatype, we have to convert it to datetime object type.**<jupyter_code>productDemandDf['Date'] = pd.to_datetime(productDemandDf['Date']) 
productDemandDf.dtypes<jupyter_output><empty_output><jupyter_text>**We have to convert order_demand to numeric type. So, we want to see whether all the values are convertible or not.**<jupyter_code>productDemandDf[pd.to_numeric(productDemandDf['Order_Demand'], errors='coerce').isnull()].head(5)<jupyter_output><empty_output><jupyter_text>**There are many values that are enclosed in parenthesis. So, we remove the parenthesis and convert string values to numeric.**<jupyter_code>productDemandDf['Order_Demand'] = productDemandDf['Order_Demand'].apply(lambda x: x.replace('(','').replace(')',''))
productDemandDf['Order_Demand'] = pd.to_numeric(productDemandDf['Order_Demand'], errors='coerce')
productDemandDf.isnull().sum()
nullDates = productDemandDf[productDemandDf['Date'].isnull()]
nullDates['Product_Category'].unique()<jupyter_output><empty_output><jupyter_text>**There are number of rows where date is null. Since we are doing time series analysis, we cannot perform prediction on the section with null dates. So, we remove rows with null dates.**<jupyter_code>cleanProductDemandDf = productDemandDf.dropna()
cleanProductDemandDf.head(5)<jupyter_output><empty_output><jupyter_text>**We want to see the total number of products ordered monthly and also the number of unique product codes.**<jupyter_code>monthlystats = cleanProductDemandDf.groupby(pd.Grouper(key='Date', freq="M")).agg({'Order_Demand':'sum','Product_Code':'nunique'}).reset_index()
monthlystats.head(5)
plt.scatter(monthlystats['Date'].values,monthlystats['Order_Demand'].values)
plt.xlabel('month')
plt.ylabel('aggregate order demand')
plt.show()
plt.scatter(monthlystats['Date'].values,monthlystats['Product_Code'].values)
plt.xlabel('month')
plt.ylabel('no of product codes')
plt.show()<jupyter_output><empty_output><jupyter_text>**From both plots, we see that the data for demand and unique product codes seem to be incomplete in 2011 because it is not consistent with the remaining trend. So, we discard those data. Also the data for 2017 also seem to be incomplete and we discard this too.**<jupyter_code>filteredProdDemand = productDemandDf[(productDemandDf.Date>pd.Timestamp('2011-12-31')) & \
                                     (productDemandDf.Date<pd.Timestamp('2017-1-1'))]
print(filteredProdDemand['Date'].min())
print(filteredProdDemand['Date'].max())<jupyter_output>2012-01-01 00:00:00
2016-12-30 00:00:00<jupyter_text>**Lets choose one product and apply SARIMAX to predict the demand.**<jupyter_code>prod1000 = filteredProdDemand[filteredProdDemand['Product_Code']=='Product_1000'].copy() 
prod1000['Warehouse'] = prod1000['Warehouse'].astype('category')
plt.scatter(prod1000['Date'].values, prod1000['Order_Demand'].values,c=prod1000.Warehouse.cat.codes)
plt.xlabel('weeks')
plt.ylabel('order demand for individual warehouse')
print(prod1000.shape)
print(prod1000.describe())<jupyter_output>(751, 5)
        Order_Demand
count     751.000000
mean     4252.263648
std     11688.666941
min         0.000000
25%       150.000000
50%       500.000000
75%      3000.000000
max    152150.000000<jupyter_text>**There are multiple warehouses associated with a single product and the trend of the product demand varies with each warehouse as shown in the scatter plot. So, we will use SARIMAX model for individual warehouse.**<jupyter_code>prod1000byWare = prod1000.groupby('Warehouse')['Order_Demand'].agg(median='median', mean='mean', std='std')
print(prod1000byWare)

fig,axes = plt.subplots(1,3)
for i, ware in enumerate(prod1000byWare.index.categories.values):
    axes[i].hist(prod1000[prod1000['Warehouse']==ware].Order_Demand, bins=50)
    axes[i].set_xlabel('order_Demand')
    axes[i].set_ylabel('count of particular order demand')
plt.show()<jupyter_output>median         mean           std
Warehouse                                   
Whse_A        500  4520.295983  13773.216544
Whse_J        150   228.915663    315.660407
Whse_S       2600  5314.615385   7640.039957<jupyter_text>**The statistics and the histogram show that the data distribution is variant for different warehouse.** **We have grouped the product demand of the selected product of selected warehouse by week since grouping by month gave few observation points.**<jupyter_code>prod1000byWare1 = prod1000.groupby([pd.Grouper(key='Date', freq='W'),pd.Grouper('Warehouse')])['Order_Demand'].sum().reset_index()
prod1000byWare1.head(5)<jupyter_output><empty_output><jupyter_text>**It is seen that the collection of data across the warehouse is not uniform since there are NaN values for some date intervals. We treat those NaN values as missing values and want our model to estimate those values.**<jupyter_code>weeklystats1000 = []
fig, ax = plt.subplots(3,1,figsize=(12,6))
for i,ware in enumerate(prod1000byWare1.Warehouse.unique()):
    value=prod1000byWare1[prod1000byWare1['Warehouse']==ware]   
    weeklystats1000.insert(i,value)
    ax[i].xaxis.grid()
    ax[i].plot(weeklystats1000[i]['Order_Demand'], 'r+')
plt.suptitle('Weekly order demand time series for 3 warehouse of product1000')
plt.show()<jupyter_output><empty_output><jupyter_text>**SARIMAX has three hyperparameters in general p,d and q, where p is the number of lag observations included in the model, also called the lag order., d is the number of times that the raw observations are differenced, also called the degree of differencing and q is the size of the moving average window, also called the order of moving average.****We have kept d =1 and found out the other 2 hyperparameters finding the lowest value of AIC(Akaike Information Criteria).**<jupyter_code>def getSARIMAXparam(data):
    
    aic = pd.DataFrame(np.zeros((6,6), dtype=float))
    warnings.simplefilter('ignore')
    # Iterate over all ARMA(p,q) models with p,q in [0,6]
    for p in range(6):
        for q in range(6):
            if p == 0 and q == 0:
                continue
                # Estimate the model with no missing datapoints
            mod = sm.tsa.statespace.SARIMAX(data, order=(p,1,q), enforce_invertibility=False)
            try:
                res = mod.fit(disp=False)
                aic.iloc[p,q] = res.aic
            except:
                aic.iloc[p,q] = np.nan
    return aic
dataRaw = []
aicList = []
for i, data in enumerate(weeklystats1000):
    dataRaw.insert(i, data['Order_Demand'].values)
    aicList.insert(i,getSARIMAXparam(dataRaw[i]))    
aicList[0]
aicList[1]
aicList[2]
hyperParam = [(2,1,5),(1,1,5),(1,1,1)]<jupyter_output><empty_output><jupyter_text>**Now,we selected hyperparameters looking at the tables above and we are going to apply SARIMAX to data of all 3 warehouse.**<jupyter_code>model = []
result = []
for i,param in enumerate(hyperParam):
    model.insert(i,sm.tsa.statespace.SARIMAX(dataRaw[i], order=param ,enforce_invertibility=False))
    result.insert(i,model[i].fit(disp=False))
predict = []
idx = []
predict_ci = []
fig, ax = plt.subplots(3, 1, figsize=(12,6))
for i in range(len(model)):
    nforecast = 20
    predict.insert(i, result[i].get_prediction(end=model[i].nobs + nforecast))
    idx.insert(i, np.arange(len(predict[i].predicted_mean)))
    predict_ci.insert(i, predict[i].conf_int(alpha=0.5))
    ax[i].xaxis.grid()
    ax[i].plot(dataRaw[i], 'r+')
    ax[i].plot(idx[i][:-nforecast], predict[i].predicted_mean[:-nforecast], 'black')
    ax[i].plot(idx[i][-nforecast:], predict[i].predicted_mean[-nforecast:], 'k--', linestyle='--', linewidth=2)
    ax[i].fill_between(idx[i], predict_ci[i][:, 0], predict_ci[i][:, 1], alpha=0.15)
plt.suptitle('Estimated timeseries using SARIMAX model')  
plt.show()<jupyter_output><empty_output><jupyter_text>**The SARIMAX model we used considers the past values to predict the future. The non stationary data is made stationary using d=1. Since, there were many missing (NaN) values in the data, this model has tried to estimate those missing values seeing the recent past value. In the graph, there are some spaces without red '+' sign. Those spaces represent the missing values. As seen in the above graph, the values are reasonably fitted by the model, if the data trend is almost stationary. However, the future values predicted are stationary around the mean.**<jupyter_code>size = int(len(dataRaw[1]) * 0.75)
train, test = dataRaw[1][0:size], dataRaw[1][size:len(dataRaw[1])]
history = [x for x in train]
predictions = list()
for t in range(len(test)):
    model = sm.tsa.statespace.SARIMAX(history, order=(1,1,5))
    model_fit = model.fit(disp=0)
    output = model_fit.forecast()
    yhat = output[0]
    predictions.append(yhat)
    obs = test[t]
    history.append(obs)
    #print('predicted=%f, expected=%f' % (yhat, obs))

notNan =  ~np.isnan(test)
error = mean_squared_error(test[notNan], np.array(predictions)[notNan])
print('Test MSE: %.3f' % error)
# plot
plt.plot(test[notNan],label='true')
plt.plot(np.array(predictions)[notNan], color='red',label = 'predicted')
plt.legend()
plt.title('true vs predicted value for next timestamp')
plt.show()<jupyter_output>Test MSE: 8013219857.711<jupyter_text>**In the above graph, although the mean square error is very high, the predicted value has tried to follow the true value. If the data distribution is stationary, this model seems work well. The model suffers, if the data has fluctuation.**  **Prophet****Now we will use prophet for all 3 cases above.**<jupyter_code>forecast = []
future = []
df = []
prophetModel = []
fig, ax = plt.subplots(3, 1, figsize=(12,6))
for i in range(len(weeklystats1000)):
    df.insert(i,weeklystats1000[i][['Date','Order_Demand']])  
    df[i].columns = ['ds','y']
    prophetModel.insert(i,Prophet(weekly_seasonality=True, daily_seasonality=False))    
    prophetModel[i].fit(df[i])
    future.insert(i, prophetModel[i].make_future_dataframe(periods=20,freq='W'))
    forecast.insert(i, prophetModel[i].predict(future[i]))
    forecast[i][['ds', 'yhat', 'yhat_lower', 'yhat_upper']]
    prophetModel[i].plot(forecast[i],ax=ax[i])
plt.suptitle('Estimated timeseries using fbprophet model')
plt.show()<jupyter_output><empty_output><jupyter_text>**As seen from the above graphs, prophet is seen to maintain the periodocity of the trend. In comparison to SARIMAX, prophet predicts the future data and fills the missing values maintaining the periodocity whereas SARIMAX predicts the future data around the mean and fills the missing values looking at the recent past trend.**<jupyter_code>size = int(len(df[1]) * 0.75)
train = df[1][0:size].reset_index()
test = df[1][size:].reset_index()
predictions = list()
for t in range(len(test)):
    model = Prophet(weekly_seasonality=True,daily_seasonality=False)
    model.fit(train)
    future = model.make_future_dataframe(periods=1,freq='W')
    output = model.predict(future)
    yhat = output['yhat'].values[0]
    predictions.append(yhat)
    ds = test['ds'].iloc[t]
    #print(type(ds))
    obs = test['y'].iloc[t]
    #print(obs)
    train.loc[len(train)] = test.iloc[t].values
   # print('predicted=%f, expected=%f' % (yhat, obs))

notNan =  ~np.isnan(test['y'].values)
error = mean_squared_error(test['y'].values[notNan], np.array(predictions)[notNan])
print('Test MSE: %.3f' % error)
# plot
plt.plot(test['y'].values[notNan],label='true')
plt.plot(np.array(predictions)[notNan], color='red',label = 'predicted')
plt.title('True vs predicted value for next timestamp')
plt.legend()
plt.show()<jupyter_output>Test MSE: 9916877484.232<jupyter_text>**We want to see how prophet predicts the value for the next timestamp. The predicted values are affected by outliers in prophet whereas in SARIMAX, prediction is quicky recovered from the outliers.**<jupyter_code>train = df[1][0:size].reset_index()
model = Prophet(weekly_seasonality=True, daily_seasonality=False)
model.fit(train)
future = model.make_future_dataframe(periods=len(test),freq='W')
output = model.predict(future)
yhat = output['yhat'].values[-len(test):]
notNan =  ~np.isnan(test['y'].values)
error = mean_squared_error(test['y'].values[notNan], yhat[notNan])
print('Test MSE: %.3f' % error)
# plot
plt.plot(test['y'].values[notNan],label='true')
plt.plot(yhat[notNan], color='red',label = 'predicted')
plt.title('True vs predicted value for future timestamps')
plt.legend()
plt.show()<jupyter_output>Test MSE: 4913383141.147
{'fraction_non_alphanumeric': 0.12493760251016188, 'fraction_numerical': 0.032803251800613276, 'mean_word_length': 5.498609823911028, 'pattern_counts': {'":': 0, '<': 64, '<?xml version=': 0, '>': 64, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '241468', 'license': 'MIT', 'n_tokens_mistral': 4921, 'n_tokens_neox': 4496, 'n_words': 1218, 'path': 'notebook/DemandPredict.ipynb', 'repo_name': 'kamala421/EECS731_WorldWideProducts', 'size': '14023'}