<jupyter_start><jupyter_text>Sales Forecasting using Walmart data set. Sales forecast is an indicator of a product quantity that is likely to be sold in a time period in a specified market at given price. In other words, accurate sales forecasting enables a company to produce the required quantity at the right time in order to meet the future demand. [1]In this project, I used a machine learning approach to forecast sales for each departemnt of each of 45 Walmart stores based on [Walmart data set](https://github.com/apadhye90/data-mining/tree/master/Kaggle%20Dataset) [2]> [1] [Sales Forecasting: Meaning, Importance and Methods](https://www.businessmanagementideas.com/sales/forecasting-sales/sales-forecasting-meaning-importance-and-methods/7122)[2] [Walmart Recruiting - Store Sales Forecasting: Overview](https://www.kaggle.com/c/walmart-recruiting-store-sales-forecasting)  Table of Contents  General Setup___<jupyter_code># Dataframes.
import pandas as pd

# Numerical arrays.
import numpy as np

# Stationarity
from statsmodels.tsa.stattools import adfuller

# Predictions
from pmdarima.arima import auto_arima
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
from sklearn.metrics import mean_squared_error
from pmdarima.arima.utils import ndiffs
from pmdarima.utils import diff_inv
import warnings
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.seasonal import STL

# Plotting.
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import seaborn as sns

%matplotlib inline
# Change style and size of plots
plt.style.use("ggplot")
plt.rcParams["figure.figsize"] = (12, 8)
plt.rcParams["figure.titlesize"] = 15<jupyter_output><empty_output><jupyter_text>Exploratory Data Analysis.___The data set consists of 4 csv files: stores, train, test, features. [3] First, I will analyse them separately. [4]> [3] [Walmart Recruiting - Store Sales Forecasting: Data Description](https://www.kaggle.com/c/walmart-recruiting-store-sales-forecasting/data)[4] [Exploratory Data Analysis(EDA): Python](https://towardsdatascience.com/exploratory-data-analysis-eda-python-87178e35b14) Stores._1. Load the file._<jupyter_code># Load the stores.csv without an index.
stores = pd.read_csv("Walmart Data Set/stores.csv", header=0)

# Display the dafaframe
stores.head()
# Display the dafaframe
stores.tail()
stores.describe()<jupyter_output><empty_output><jupyter_text>The stores file consist of information about 45 stores, including the type and size of each. We can observe that there are mainly empty values in the columns Store A, Store B and Store C.  _2. Data cleaning._The first step to cleansing the data is by checking it for empty values.<jupyter_code># Check for empty values.
stores.isnull().sum()<jupyter_output><empty_output><jupyter_text>The above confirms that only the last three columns have the empty values. Since they do not provide enough information that could be valuable in the sales forecasting, they will be removed.<jupyter_code># Remove columns with empty values.
cleaned_stores = stores.drop(['Store A','Store B','Store C'], axis=1)

# Check for empty values again.
cleaned_stores.isnull().sum()
# Calculate maximum values of each store type and set as pie sizes
sizes = cleaned_stores.groupby('Type').max().Size.values
labels = cleaned_stores.groupby('Type').max().index

# Create a figure and axis, set a title.
fig, ax = plt.subplots(figsize=(8,8))
fig.suptitle("Store Size per Type")

# Build a pie plot.
ax.pie(sizes, labels=labels, autopct='%1.1f%%', shadow=True, startangle=90)
# Equal aspect ratio ensures that pie is drawn as a circle.
ax.axis('equal') 

# Show the plot.
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Train.<jupyter_code># Load the train.csv without an index.
train = pd.read_csv("Walmart Data Set/train.csv", header=0, parse_dates=True, index_col="Date")

# Display the dafaframe.
train.head()
# Display the dafaframe.
train.tail()
train.describe()<jupyter_output><empty_output><jupyter_text>The above dataframe contains weekly sales for 45 stores per department for the period from 2010-02-05 to 2012-10-26. We can see that altogether there are 99 departments. It also includes the department number and whether the week is a special holiday week.<jupyter_code># Check for empty values.
train.isnull().sum()<jupyter_output><empty_output><jupyter_text>The dataframe doesn't have any empty values.<jupyter_code>depts = train.groupby("Store").Dept.nunique()
print(f"Maximum Departments: {np.max(depts)}\nMinimum Departments: {np.min(depts)}")<jupyter_output>Maximum Departments: 79
Minimum Departments: 61<jupyter_text>It is clear that the stores have different distribution of departments.<jupyter_code># Plot the data by year
g = sns.catplot(x='IsHoliday', y='Weekly_Sales', data=train)
g.fig.set_size_inches(10,8)

# Show the plot.
plt.title("Weekly Sales vs Holiday")
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Features.<jupyter_code># Load the features.csv without an index.
features = pd.read_csv("Walmart Data Set/features.csv", header=0, parse_dates=True, index_col="Date")

# Display the dafaframe
features.head()
# Display the dafaframe
features.tail()<jupyter_output><empty_output><jupyter_text>The features file contains additional data related to the store, department, and regional activity for the given dates. MarkDown columns are related to promotional markdowns that Walmart is running.<jupyter_code>features.describe()
# Check for empty values.
features.isnull().sum()
features.loc['2011-11-1': '2011-11-30']<jupyter_output><ipython-input-19-7cbe5f182b62>:1: FutureWarning: Value based partial slicing on non-monotonic DatetimeIndexes with non-existing keys is deprecated and will raise a KeyError in a future Version.
  features.loc['2011-11-1': '2011-11-30']<jupyter_text>We can see that MarkDowns appear from 11.11.2011.<jupyter_code>features = features.loc['2011-11-11': '2013-07-26']
features
def date_format(ax, int):
    # Make the x axis display well.
    weeks = mdates.DayLocator(int)
    h_fmt = mdates.DateFormatter('%d-%m-%Y')

    # Tick ax axis.
    ax.xaxis.set_major_locator(weeks)
    ax.xaxis.set_major_formatter(h_fmt)
# Calculate mean of columns in features.
temp = features[features["Store"]==1].groupby(["Date"]).mean()

# Create a figure and a set of axis
fig, [[ax1,ax2],[ax3,ax4]] = plt.subplots(2,2, figsize=(15,12))

# Date Format.
fig.autofmt_xdate()

# Plot the data.
ax1.plot(temp.MarkDown1)
ax1.title.set_text("MarkDown 1")
ax2.plot(temp.MarkDown2)
ax2.title.set_text("MarkDown 2")
ax3.plot(temp.MarkDown3)
ax3.title.set_text("MarkDown 3")
ax4.plot(temp.MarkDown4)
ax4.title.set_text("MarkDown 4")

# Show the plots.
plt.show()
# Fill the missing values.
features["MarkDown1"].fillna(method="ffill", inplace=True)
features["MarkDown2"].fillna(method="ffill", inplace=True)
features["MarkDown3"].fillna(method="bfill", inplace=True)
features["MarkDown4"].fillna(method="ffill", inplace=True)
# Create a figure and a set of axis
fig, [ax1,ax2] = plt.subplots(1,2)

# Date Format.
fig.autofmt_xdate()

# Plot the data.
ax1.plot(temp.CPI)
ax1.title.set_text("CPI")
ax2.plot(temp.Unemployment)
ax1.title.set_text("Unemployment")

# Show the plots.
plt.show()
# Filling missing values for CPI and unemployment
features["CPI"].fillna(method="ffill", inplace=True)
features["Unemployment"].fillna(method="ffill", inplace=True)
# Checking for missing values again.
features.isnull().sum()
# Separate test features from train.
features_train = features.loc['2011-11-11': '2012-10-26'].reset_index()
features_test = features.loc['2012-11-02': '2013-07-26'].reset_index()
features_train.head()<jupyter_output><empty_output><jupyter_text>Mergining data sets.___<jupyter_code># Merge the dataframes into one.
df = train.merge(features_train, on=['Store', 'IsHoliday', 'Date'],how='left').dropna()
df = df.merge(cleaned_stores, on=['Store'], how='left')

# Timedelta
#df['Date'] = df['Date'] + pd.to_timedelta(df.groupby('Date').cumcount(), unit='m')
df=df.set_index(['Date'])

# Display the new dataframe with a date index.
df
df.index.duplicated().sum()
df.isnull().sum()
# Create a figure and axis.
fig, ax = plt.subplots(figsize=(13,10))
fig.suptitle("Weekly Sales")

# Date formatter.
date_format(ax,10)
fig.autofmt_xdate()

# Plot the sales
sns.scatterplot(data=df, x=df.index, y="Weekly_Sales", ax=ax)

# Show the plot.
plt.tight_layout()
plt.show()
# Create a figure and axis.
fig, ax = plt.subplots(figsize=(13,10))
fig.suptitle("Correlation Map")

# Plot the correlation
corr = df.corr()
cmap = sns.diverging_palette(100, 275, as_cmap=True)
sns.heatmap(corr, cmap=cmap, center=0, annot=True,square=True,cbar_kws={"shrink": .7}, ax=ax)

# Set a title and rotate x axis 
plt.xticks(rotation=45)
plt.tight_layout()

# Show the plot.
plt.show()<jupyter_output><empty_output><jupyter_text>The heatmap shows that weekly sales are dependent on:* Store Size* Department* MarkDown5, MarkDown1, MarkDown3, MarkDown4, MarkDown2* IsHoliday      Holidays & Markdowns<jupyter_code># Create a figure and axis.
fig, ax = plt.subplots(figsize=(20,8))
fig.suptitle("Weekly Sales per Store")

# Plot the data by year
per_store = df.groupby("Store").agg({"Weekly_Sales": "sum"})
per_store.sort_values("Weekly_Sales").plot.bar(ax=ax)

# Set a title and change rotation of xticks back to 0 
plt.title("2011-2012")
plt.xticks(rotation=0)

# Show the plot.
plt.show()
"""
def test_stationarity(ts):
    # Create a figure and axis.
    fig, ax = plt.subplots()
    fig.suptitle('Rolling Mean & Standard Deviation')

    # Date Format.
    date_format(ax,10)
    fig.autofmt_xdate()

    # Determing rolling statistics
    rolmean = ts["Weekly_Sales"].rolling(12).mean()
    rolstd = ts["Weekly_Sales"].rolling(12).std()

    # Plot rolling statistics:
    sns.lineplot(data=ts, x="Date", y="Weekly_Sales", ax=ax, label='Original')
    sns.lineplot(data=ts, x="Date", y=rolmean, ax=ax, label='Rolling Mean')
    sns.lineplot(data=ts, x="Date", y=rolstd, ax=ax, label = 'Rolling Std')

    # Show the plot
    plt.show()

    #Perform Dickey-Fuller test:
    print('Results of Dickey-Fuller Test:')
    dftest = adfuller(ts, autolag='AIC')
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
    for key,value in dftest[4].items():
        dfoutput[f'Critical Value ({key})'] = value
    print(dfoutput)

sales_df = pd.DataFrame(df["Weekly_Sales"])
test_stationarity(sales_df)
"""<jupyter_output><empty_output><jupyter_text>Preprocessing<jupyter_code># Change to pivot table per department and tore
sales = pd.pivot_table(df,index=df.index,columns=[df.Store, df.Dept], values='Weekly_Sales',aggfunc=np.sum)

# Stack the pivot table and remove rows with more than 10% (apprx 5) missing values
sales = sales.dropna(thresh=len(sales) - 5, axis=0)
sales = sales.fillna(method="bfill")
sales = sales.fillna(method="ffill")
sales.isnull().sum().any()
sales.head()
res = STL(sales.values.flatten(), period=7).fit()
res.plot()
plt.show()<jupyter_output><empty_output><jupyter_text>ARIMA model.___<jupyter_code>def get_dept():
    depts_sales = []
    for i,j in sales:
        dept_sales = sales.loc[:, (i, j)]
        depts_sales.append(dept_sales)
    return depts_sales
def split(dept_sales):
    size = int(len(dept_sales) * 0.7)
    train, val = dept_sales[0:size], dept_sales[size:]
    return train, val
def adf(dept_sales):
    dftest = adfuller(dept_sales, autolag='AIC')
    values = []
    for key,value in dftest[4].items():
        values.append(value)
    return dftest[0], values[0]
def make_stationary(dept_sales):
    dept_sales_diff = dept_sales.diff().dropna()
    return dept_sales_diff
def inv_diff (df_orig_column,df_diff_column, periods):
    # Generate np.array for the diff_inv function - it includes first n values(n = 
    # periods) of original data & further diff values of given periods
    value = np.array(df_orig_column[:periods].tolist()+df_diff_column[periods:].tolist())

    # Generate np.array with inverse diff
    inv_diff_vals = diff_inv(value, periods,1 )[periods:]
    return inv_diff_vals
def arima(dept_sales, train, val):
    # Create a model
    model = auto_arima(dept_sales, start_p=0, start_q=0, start_P=0, start_Q=0,
                        test='adf', # use Augmented Dickey-Fuller(adftest) test to find optimal 'd'
                      trend="ct", stationary=True,
                      stepwise=False)
    model.fit(train)
    preds = model.predict(n_periods=len(val))
    val_vs_forecast = plot(preds, val)
    return val_vs_forecast
def plot(preds, val):
    forecast_df = pd.DataFrame(preds,index = val.index,columns=['Prediction'])
    
    val_vs_forecast = pd.concat([val,forecast_df],axis=1)
    val_vs_forecast.plot()
    plt.title(f'RMSE: {rmse(val, forecast_df):.2f}')    
    return val_vs_forecast
def rmse(val, preds):
    rmse = np.sqrt(mean_squared_error(val, preds))
    return rmse
def predictions_per_dept(dept_sales):
    train, val = split(dept_sales)
    #test_stat, crit_val = adf(dept_sales)
    #if test_stat > crit_val:
    #    dept_sales_diff = make_stationary(dept_sales)
    #    train_diff, val_diff = split(dept_sales_diff)
     #   val_vs_forecast = arima(dept_sales_diff, train_diff, val_diff)
    val_vs_forecast = arima(dept_sales, train, val)
    return val_vs_forecast
predictions = pd.DataFrame()
depts_sales = get_dept()
predictions = predictions.merge(predictions_per_dept(depts_sales[9]), on="Date", how="left")
predictions.T
# Future forecast
predictions = pd.DataFrame()
warnings.filterwarnings("ignore")
depts_sales = get_dept()
for i in range(0, 3):
    preds = predictions_per_dept(depts_sales[i])
    predictions = predictions.append(preds)
predictions
# Create a figure and axis.
fig, ax = plt.subplots()

ax.plot(val, label='Original')
ax.plot(preds, label='Predicted', color='b')

plt.legend()
plt.show()
train, val = split(pd.DataFrame(depts_sales[9]))
# evaluate an ARIMA model for a given order (p,d,q)
def evaluate_arima_model(order):
    # prepare training dataset
    model = ARIMA(train, order=order).fit()
    predictions = model.predict(start=1, end=len(val), typ='levels')
    # calculate out of sample error
    rmse = np.sqrt(mean_squared_error(val, predictions))
    return rmse
# evaluate combinations of p, d and q values for an ARIMA model
def evaluate_models(p_values, d_values, q_values):
    best_score, best_cfg = float("inf"), None
    for p in p_values:
        for d in d_values:
            for q in q_values:
                order = (p,d,q)
                rmse = evaluate_arima_model(order)
                if rmse < best_score:
                    best_score, best_cfg = rmse, order
    print(f'Best ARIMA {best_cfg}: RMSE={best_score:.3f}')

# evaluate parameters
p_values = [0, 1, 2, 4, 6, 8, 10]
d_values = range(0, 3)
q_values = range(0, 3)
warnings.filterwarnings("ignore")
evaluate_models(p_values, d_values, q_values)<jupyter_output><empty_output><jupyter_text>Test data.<jupyter_code># Load the train.csv without an index.
test_data = pd.read_csv("Walmart Data Set/test.csv", header=0, parse_dates=True, index_col="Date")

# Display the dafaframe.
test_data.head()
test_data.tail()
test_df = test_data.merge(features_test, on=['Store', 'IsHoliday', 'Date'],how='left')
test_df.index = test_df.Date
test_df = test_df.drop("Date", axis=1)
test_df.head()
test_df.tail()
test_df.isnull().sum()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12726222564497497, 'fraction_numerical': 0.01893210114234373, 'mean_word_length': 5.196023856858847, 'pattern_counts': {'":': 2, '<': 65, '<?xml version=': 0, '>': 67, 'https://': 5, 'lorem ipsum': 0, 'www.': 3, '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': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '273337', 'license': 'MIT', 'n_tokens_mistral': 5686, 'n_tokens_neox': 5352, 'n_words': 1347, 'path': 'Sales Forecasting using Walmart Dataset.ipynb', 'repo_name': 'olgarozhdestvina/Sales-Forecasting', 'size': '15582'}