<jupyter_start><jupyter_text>The dataset contains historical product demand for a manufacturing company with footprints globally. The company provides thousands of products within dozens of product categories. There are four central warehouses to ship products within the region it is responsible for. Since the products are manufactured in different locations all over the world, it normally takes more than one month to ship products via ocean to different central warehouses. If forecasts for each product in different central with reasonable accuracy for the monthly demand for month after next can be achieved, it would be beneficial to the company in multiple ways.The information contained in the dataset was encoded to protect the company's identity. Date variable must be formated as datetime for time series analysis.<jupyter_code>Datos_datetime=Datos
Datos_datetime["Date"]=pd.to_datetime(Datos["Date"],format="%Y/%m/%d")
Datos_datetime.dtypes<jupyter_output><empty_output><jupyter_text>Checking for missing values<jupyter_code>print(Datos_datetime.isnull().sum())
print ('Null to Dataset Ratio for "Dates" Column '': ',Datos_datetime.isnull().sum()[3]/Datos_datetime.shape[0]*100)<jupyter_output>Product_Code            0
Warehouse               0
Product_Category        0
Date                11239
Order_Demand            0
dtype: int64
Null to Dataset Ratio for "Dates" Column :  1.0718355863910547<jupyter_text>The null data will be removed due to it being around 1% of the data.<jupyter_code>Datos_datetime.dropna(axis=0, inplace=True)
Datos_datetime.reset_index(drop=True)
Datos_datetime.sort_values('Date')[1:50]
Datos_datetime['Order_Demand']=Datos_datetime['Order_Demand'].str.replace('(',"")
Datos_datetime['Order_Demand']=Datos_datetime['Order_Demand'].str.replace(')',"")
Datos_clean=Datos_datetime
Datos_clean.head()
Datos_clean['Order_Demand'] = Datos_clean['Order_Demand'].astype('int64')
Datos_clean.dtypes
print(Datos_clean.isnull().sum())
print ('Null to Dataset Ratio for "Dates" Column '': ',Datos_clean.isnull().sum()[3]/Datos_clean.shape[0]*100)
Datos_clean['Date'].min() , Datos_clean['Date'].max()<jupyter_output><empty_output><jupyter_text>The time series goes from January 2011 to January 2017.<jupyter_code>Datos_clean['Warehouse'].value_counts().sort_values(ascending=False)
Datos_WH=Datos_clean.groupby('Warehouse').sum().sort_values('Order_Demand', ascending = False)
Datos_WH
rcParams['figure.figsize']=10,10 
logDataWH=sn.boxplot(Datos_clean['Warehouse'],np.log1p(Datos_clean['Order_Demand'])) 
del  logDataWH<jupyter_output>E:\Anaconda\lib\site-packages\seaborn\_decorators.py:36: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.
  warnings.warn(<jupyter_text>Product Code exploration<jupyter_code>print(len(Datos_clean['Product_Code'].value_counts()))<jupyter_output>2160<jupyter_text>There are 2160 different products Product Category exploration<jupyter_code>print("Product Categories:")
print(len(Datos_clean['Product_Category'].value_counts()))
rcParams['figure.figsize']=50,14
sn.countplot(Datos_clean['Product_Category'].sort_values(ascending=True))<jupyter_output>Product Categories:
33<jupyter_text>There are 33 existing product categories.<jupyter_code>rcParams['figure.figsize'] = 50,12
df_temp = Datos_clean.sample(n=20000).reset_index()
fig5 = sn.boxplot( df_temp['Product_Category'].sort_values(),np.log1p(df_temp['Order_Demand']))
del df_temp, fig5<jupyter_output>E:\Anaconda\lib\site-packages\seaborn\_decorators.py:36: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.
  warnings.warn(<jupyter_text>Due to the existance of too many distinct products and the differences in the product category means, the demand forecasting will be done according to product categories.<jupyter_code>Datos_PC=Datos_clean.groupby('Product_Category').count()
Datos_PC=Datos_PC.rename(columns={"Warehouse":"Frequency"})
Datos_PC.drop(columns={"Product_Code","Date","Order_Demand"})
DatosPC_sorted=Datos_PC["Frequency"].sort_values(ascending=False)
DatosPC_sorted<jupyter_output><empty_output><jupyter_text>The categories of choice for the predictions done in this study were 019, 005 and 001. Since these three are the most frequent categories present on this dataset. Here we separate the dataset in 3 depending on the product category.<jupyter_code>Datos_Cat019=Datos_clean[Datos_clean["Product_Category"]=="Category_019"]
Datos_Cat019=Datos_Cat019[["Date", "Order_Demand"]]
Datos_Cat001=Datos_clean[Datos_clean["Product_Category"]=="Category_001"]
Datos_Cat001=Datos_Cat001[["Date", "Order_Demand"]]
Datos_Cat005=Datos_clean[Datos_clean["Product_Category"]=="Category_005"]
Datos_Cat005=Datos_Cat005[["Date", "Order_Demand"]]
TS_019=Datos_Cat019.set_index("Date")
TS_005=Datos_Cat005.set_index("Date")
TS_001=Datos_Cat001.set_index("Date")
TS_019 = TS_019['Order_Demand'].resample('MS').mean()
TS_019 = TS_019.fillna(TS_019.bfill())
TS_005 = TS_005['Order_Demand'].resample('MS').mean()
TS_005 = TS_005.fillna(TS_005.bfill())
TS_001 = TS_001['Order_Demand'].resample('MS').mean()
TS_001 = TS_001.fillna(TS_001.bfill())<jupyter_output><empty_output><jupyter_text>Here we separate the time series in training and testing. Having 70% of the data for training and the 30% for testing.<jupyter_code>lim_019=round(len(TS_019)*0.7)
TS_019_Training=TS_019[1:lim_019]
TS_019_Testing=TS_019[lim_019+1:]

lim_005=round(len(TS_005)*0.7)
TS_005_Training=TS_005[1:lim_005]
TS_005_Testing=TS_005[lim_005+1:]

lim_001=round(len(TS_005)*0.7)
TS_001_Training=TS_001[1:lim_001]
TS_001_Testing=TS_001[lim_001+1:]

#Defining function for Time series decomposition
def decompose(data):
    rcParams['figure.figsize'] = 20, 10
    decomposition = sm.tsa.seasonal_decompose(data, model='additive')
    fig = decomposition.plot()
    return plt.show()<jupyter_output><empty_output><jupyter_text>Category 019 Time Series Decomposition<jupyter_code>decompose(TS_019_Training)<jupyter_output><empty_output><jupyter_text>Category 005 Time Series Decomposition<jupyter_code>decompose(TS_005_Training)<jupyter_output><empty_output><jupyter_text>Category 001 Time Series Decomposition<jupyter_code>decompose(TS_001_Training)<jupyter_output><empty_output><jupyter_text>Dickey Fuller stationarity test<jupyter_code>#Defining function for Dickey Fuller test
def DF_test(data):
    print("Results of Dickey Fuller Test:")
    dftest=adfuller(data,autolag="AIC")
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
    print(dfoutput)
DF_test(TS_019_Training)
DF_test(TS_005_Training)
DF_test(TS_001_Training)<jupyter_output>Category 019 Time Series
Results of Dickey Fuller Test:
Test Statistic                 -5.103645
p-value                         0.000014
#Lags Used                      0.000000
Number of Observations Used    46.000000
dtype: float64<jupyter_text>ARIMA models Auto-ARIMA<jupyter_code>#Defining function for Auto Arima modeling    
def AutoArima(data):
    return auto_arima(data, start_p=1, start_q=1,
                          max_p=5, max_q=5,
                          start_P=0, seasonal=False,
                          d=None, D=1, trace=True,
                          error_action='ignore',   
                          suppress_warnings=True,  
                          stepwise=True)           
model_019_autoarima=AutoArima(TS_019_Training)
model_019_autoarima.summary()
model_005_autoarima=AutoArima(TS_005_Training)
model_005_autoarima.summary()
model_001_autoarima=AutoArima(TS_001_Training)
model_001_autoarima.summary()<jupyter_output>Performing stepwise search to minimize aic
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.10 sec
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=162.075, Time=0.01 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=158.476, Time=0.01 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.05 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=160.328, Time=0.01 sec
 ARIMA(2,1,0)(0,0,0)[0] intercept   : AIC=158.631, Time=0.02 sec
 ARIMA(2,1,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.10 sec
 ARIMA(1,1,0)(0,0,0)[0]             : AIC=157.029, Time=0.01 sec
 ARIMA(2,1,0)(0,0,0)[0]             : AIC=157.461, Time=0.01 sec
 ARIMA(1,1,1)(0,0,0)[0]             : AIC=153.612, Time=0.02 sec
 ARIMA(0,1,1)(0,0,0)[0]             : AIC=153.240, Time=0.01 sec
 ARIMA(0,1,2)(0,0,0)[0]             : AIC=153.710, Time=0.02 sec
 ARIMA(1,1,2)(0,0,0)[0]             : AIC=inf, Time=0.08 sec

Best model:  ARIMA(0,1,1)(0,0,0)[0]          
Total fit time: 0.446 seconds<jupyter_text>ARIMA Models Visualization<jupyter_code>Model_019_arima=ARIMA(TS_019_Training,order=model_019_autoarima.order)
Results_M_019=Model_019_arima.fit()
plt.plot(TS_019_Training)
plt.plot(Results_M_019.fittedvalues, color="red")
plt.title("Modelo ARMA Cat 019")
Model_005_arima=ARIMA(TS_005_Training,order=model_005_autoarima.order)
Results_M_005=Model_005_arima.fit()
plt.plot(TS_005_Training)
plt.plot(Results_M_005.fittedvalues[1:], color="red")
plt.title("Modelo ARMA Cat 005")
Model_001_arima=ARIMA(TS_001_Training,order=model_001_autoarima.order)
Results_M_001=Model_001_arima.fit()
plt.plot(TS_001_Training)
plt.plot(Results_M_001.fittedvalues[1:], color="red")
plt.title("Modelo ARMA Cat 001")<jupyter_output><empty_output><jupyter_text>ARIMA Models Forecats<jupyter_code>Forecast_Arima_019=model_019_autoarima.predict(n_periods=len(TS_019_Testing))
Forecast_Arima_019
Forecast_Arima_005=model_005_autoarima.predict(n_periods=len(TS_005_Testing))
Forecast_Arima_005
Forecast_Arima_001=model_001_autoarima.predict(n_periods=len(TS_001_Testing))
Forecast_Arima_001<jupyter_output><empty_output><jupyter_text>Error metrics for ARIMA models<jupyter_code>MAE_ARIMA_019=mean_absolute_error(Forecast_Arima_019,TS_019_Testing)
MAE_ARIMA_005=mean_absolute_error(Forecast_Arima_005,TS_005_Testing)
MAE_ARIMA_001=mean_absolute_error(Forecast_Arima_001,TS_001_Testing)
print("MAE ARIMA Model Category 019: ",MAE_ARIMA_019)
print("MAE ARIMA Model Category 005: ",MAE_ARIMA_005)
print("MAE ARIMA Model Category 001: ",MAE_ARIMA_001)<jupyter_output>MAE ARIMA Model Category 019:  635.6685190444092
MAE ARIMA Model Category 005:  224.92575675655718
MAE ARIMA Model Category 001:  2.076196988321565<jupyter_text>XGBoost Models<jupyter_code>#Defining function to transform Time Series object into DataFrame
def TS_to_dataframe(time_series):
    df=pd.DataFrame(time_series)
    df.reset_index(inplace=True)
    return df

TS_019_Training_df=TS_to_dataframe(TS_019_Training)
TS_005_Training_df=TS_to_dataframe(TS_005_Training)
TS_001_Training_df=TS_to_dataframe(TS_001_Training)

TS_019_Testing_df=TS_to_dataframe(TS_019_Testing)
TS_005_Testing_df=TS_to_dataframe(TS_005_Testing)
TS_001_Testing_df=TS_to_dataframe(TS_001_Testing)
def create_features(df, label=None):
    """
    Creates time series features from Datetime index
    """
    df['hour'] = df['Date'].dt.hour
    df['dayofweek'] = df['Date'].dt.dayofweek
    df['quarter'] = df['Date'].dt.quarter
    df['month'] = df['Date'].dt.month
    df['year'] = df['Date'].dt.year
    df['dayofyear'] = df['Date'].dt.dayofyear
    df['dayofmonth'] = df['Date'].dt.day
    df['weekofyear'] = df['Date'].dt.weekofyear
    
    X = df[['hour','dayofweek','quarter','month','year',
           'dayofyear','dayofmonth','weekofyear']]
    if label:
        y = df[label]
        return X, y
    return X
X_Train_019,Y_Train_019=create_features(TS_019_Training_df,"Order_Demand")
X_Test_019,Y_Test_019=create_features(TS_019_Testing_df,"Order_Demand")

X_Train_005,Y_Train_005=create_features(TS_005_Training_df,"Order_Demand")
X_Test_005,Y_Test_005=create_features(TS_005_Testing_df,"Order_Demand")

X_Train_001,Y_Train_001=create_features(TS_001_Training_df,"Order_Demand")
X_Test_001,Y_Test_001=create_features(TS_001_Testing_df,"Order_Demand")

XGB_019 = XGBRegressor(n_estimators=500)
XGB_019.fit(X_Train_019,Y_Train_019,eval_metric="mae")
Forecast_XGB_019=XGB_019.predict(X_Test_019)
print(Forecast_XGB_019)

plot_importance(XGB_019, height=0.9)
XGB_005 = XGBRegressor(n_estimators=500)
XGB_005.fit(X_Train_005,Y_Train_005,eval_metric="mae")
Forecast_XGB_005=XGB_005.predict(X_Test_005)
print(Forecast_XGB_005)
plot_importance(XGB_005, height=0.9)
XGB_001 = XGBRegressor(n_estimators=500)
XGB_001.fit(X_Train_001,Y_Train_001,eval_metric="mae")
Forecast_XGB_001=XGB_001.predict(X_Test_001)
print(Forecast_XGB_001)
plot_importance(XGB_001, height=0.9)
MAE_XGB_019=mean_absolute_error(Forecast_XGB_019,Y_Test_019)
MAE_XGB_005=mean_absolute_error(Forecast_XGB_005,Y_Test_005)
MAE_XGB_001=mean_absolute_error(Forecast_XGB_001,Y_Test_001)
print("MAE XGBoost Model Category 019: ",MAE_XGB_019)
print("MAE XGBoost Model Category 005: ",MAE_XGB_005)
print("MAE XGBoost Model Category 001: ",MAE_XGB_001)<jupyter_output>MAE XGBoost Model Category 019:  633.9247722517326
MAE XGBoost Model Category 005:  266.29841274404083
MAE XGBoost Model Category 001:  2.2267223921128503<jupyter_text>Model Ensemble The method for combining both models into a single model ensemble is averaging the predictions of both models.<jupyter_code>ensemble_prediction_019=(Forecast_XGB_019+Forecast_Arima_019)/2
ensemble_prediction_005=(Forecast_XGB_005+Forecast_Arima_005)/2
ensemble_prediction_001=(Forecast_XGB_001+Forecast_Arima_001)/2<jupyter_output><empty_output><jupyter_text>Error metrics for model ensemble<jupyter_code>MAE_Ensemble_019=mean_absolute_error(ensemble_prediction_019,Y_Test_019)
MAE_Ensemble_005=mean_absolute_error(ensemble_prediction_005,Y_Test_005)
MAE_Ensemble_001=mean_absolute_error(ensemble_prediction_001,Y_Test_001)
print("MAE Model Ensemble Category 019: ",MAE_Ensemble_019)
print("MAE Model Ensemble Category 005: ",MAE_Ensemble_005)
print("MAE Model Ensemble Category 001: ",MAE_Ensemble_001)<jupyter_output>MAE Model Ensemble Category 019:  561.6406747496818
MAE Model Ensemble Category 005:  238.79895122578375
MAE Model Ensemble Category 001:  2.1419957904003644
{'fraction_non_alphanumeric': 0.11345831001678791, 'fraction_numerical': 0.08568830442081701, 'mean_word_length': 7.346176298890835, 'pattern_counts': {'":': 1, '<': 75, '<?xml version=': 0, '>': 75, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '389404', 'license': 'MIT', 'n_tokens_mistral': 6071, 'n_tokens_neox': 5269, 'n_words': 875, 'path': 'ARIMA+XGBoost.ipynb', 'repo_name': 'nethKen1/Gradient-Boosted-ARIMA-For-Time-Series-Forecasting', 'size': '14296'}