Defining the Business Problem¶

HappyMoonVC is a Venture Capital (VC) firm.

Venture Capital is a class of investment fund focused in growth capital for medium sized companies which already possess a portfolio of customers as well as reveneu, but which still need to gain scale. VCs invest in Startups and young companies with high growth potential.

HappyMoonVC has already invested in a new promissing Startup called GoInRail.

GoInRail is an innovative transportation method which utilized jet propulsion technology to launch a drone-lie vehicle over rails (which maintain the drone charged with electricity), and thus transport people in high speed!

Now HappyMoonVC is analyzing an complementary business, an e-commerce platform capable of supplying goods for GoInRail's clients. The startup offers a touch-screen monitor inside the GoInRail which allows travelers to purchase while onboard, with the option to retrieve the goods in the station where they exit or to receive the goods at home.

(If you find this too futurist, know that it already exists in Japan's metro!)

To decide if the investment is worth HappyMoonVC must predict the sales volumes.

Load Packages¶

In [0]:

!pip install pmdarima

Collecting pmdarima

  Downloading https://files.pythonhosted.org/packages/ff/07/7c173cc4fee44ebd62ddf03b3de84c4f151ec23facdf16baf58b8d02784c/pmdarima-1.6.0-cp36-cp36m-manylinux1_x86_64.whl (1.5MB)

     |████████████████████████████████| 1.5MB 4.8MB/s

Requirement already satisfied: Cython>=0.29 in /usr/local/lib/python3.6/dist-packages (from pmdarima) (0.29.17)

Requirement already satisfied: numpy>=1.17.3 in /usr/local/lib/python3.6/dist-packages (from pmdarima) (1.18.4)

Requirement already satisfied: joblib>=0.11 in /usr/local/lib/python3.6/dist-packages (from pmdarima) (0.14.1)

Requirement already satisfied: scikit-learn>=0.22 in /usr/local/lib/python3.6/dist-packages (from pmdarima) (0.22.2.post1)

Requirement already satisfied: pandas>=0.19 in /usr/local/lib/python3.6/dist-packages (from pmdarima) (1.0.3)

Requirement already satisfied: statsmodels>=0.10.2 in /usr/local/lib/python3.6/dist-packages (from pmdarima) (0.10.2)

Requirement already satisfied: scipy>=1.3.2 in /usr/local/lib/python3.6/dist-packages (from pmdarima) (1.4.1)

Requirement already satisfied: urllib3 in /usr/local/lib/python3.6/dist-packages (from pmdarima) (1.24.3)

Requirement already satisfied: python-dateutil>=2.6.1 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.19->pmdarima) (2.8.1)

Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.19->pmdarima) (2018.9)

Requirement already satisfied: patsy>=0.4.0 in /usr/local/lib/python3.6/dist-packages (from statsmodels>=0.10.2->pmdarima) (0.5.1)

Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.6/dist-packages (from python-dateutil>=2.6.1->pandas>=0.19->pmdarima) (1.12.0)

Installing collected packages: pmdarima

Successfully installed pmdarima-1.6.0

In [0]:

# Deactivating the multiple warning messages produced by the newest versions of Pandas and Matplotlib.

import sys

import warnings

import matplotlib.cbook

if not sys.warnoptions:

    warnings.simplefilter("ignore")

warnings.simplefilter(action='ignore', category=FutureWarning)

warnings.filterwarnings("ignore", category=FutureWarning)

warnings.filterwarnings("ignore", category=matplotlib.cbook.mplDeprecation)

# Data manipulation imports

import numpy as np

import pandas as pd

import itertools

from pandas import Series

import holidays

# Data visualization imports

import matplotlib.pyplot as plt

import matplotlib as m

import seaborn as sns

import plotly as py

import plotly.express as px

import plotly.graph_objs as go

from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot

# Predictive modeling imports

import statsmodels

import statsmodels.api as sm

import statsmodels.tsa.api as smt

import statsmodels.stats as sms

import pmdarima as pm

import scipy

import scipy.stats as scs

from statsmodels.graphics import tsaplots

from statsmodels.tsa.seasonal import seasonal_decompose

from statsmodels.tsa.stattools import adfuller

from statsmodels.tsa.arima_model import ARMA

from statsmodels.tsa.arima_model import ARIMA

from statsmodels.stats.stattools import jarque_bera

import fbprophet

from fbprophet import Prophet

# Metrics and model performance imports

import math

from math import sqrt

import sklearn

from sklearn.metrics import mean_squared_error

# Graphics formatting imports

m.rcParams['axes.labelsize'] = 14

m.rcParams['xtick.labelsize'] = 12

m.rcParams['ytick.labelsize'] = 12

m.rcParams['text.color'] = 'k'

from matplotlib.pylab import rcParams

rcParams['figure.figsize'] = 20,10

matplotlib.style.use('ggplot')

%matplotlib inline

Dataset source: https://community.tableau.com/docs/DOC-1236

HappyMoonVC wants to analyze and predict all data the in the "technology" category

In [0]:

# Load data

data = pd.read_csv('https://raw.githubusercontent.com/Matheus-Schmitz/Datasets/master/tableau_superstore_sales.csv')

In [0]:

# Shape

data.shape

Out[0]:

(9994, 21)

In [0]:

# Columns

data.columns

Out[0]:

Index(['Row ID', 'Order ID', 'Order Date', 'Ship Date', 'Ship Mode',

       'Customer ID', 'Customer Name', 'Segment', 'Country', 'City', 'State',

       'Postal Code', 'Region', 'Product ID', 'Category', 'Sub-Category',

       'Product Name', 'Sales', 'Quantity', 'Discount', 'Profit'],

      dtype='object')

Exploratory Analysis¶

In [0]:

# Visualizing data

data.head()

Out[0]:

In [0]:

# Statistic summaries

data.describe()

Out[0]:

In [0]:

# Checking for missing values

data.info()

<class 'pandas.core.frame.DataFrame'>

RangeIndex: 9994 entries, 0 to 9993

Data columns (total 21 columns):

 #   Column         Non-Null Count  Dtype 

---  ------         --------------  ----- 

 0   Row ID         9994 non-null   int64 

 1   Order ID       9994 non-null   object

 2   Order Date     9994 non-null   object

 3   Ship Date      9994 non-null   object

 4   Ship Mode      9994 non-null   object

 5   Customer ID    9994 non-null   object

 6   Customer Name  9994 non-null   object

 7   Segment        9994 non-null   object

 8   Country        9994 non-null   object

 9   City           9994 non-null   object

 10  State          9994 non-null   object

 11  Postal Code    9994 non-null   int64 

 12  Region         9994 non-null   object

 13  Product ID     9994 non-null   object

 14  Category       9994 non-null   object

 15  Sub-Category   9994 non-null   object

 16  Product Name   9994 non-null   object

 17  Sales          9994 non-null   float64

 18  Quantity       9994 non-null   int64 

 19  Discount       9994 non-null   float64

 20  Profit         9994 non-null   float64

dtypes: float64(3), int64(3), object(15)

memory usage: 1.6+ MB

In [0]:

# Setting column names to lowercase

data.columns = map(str.lower, data.columns)

In [0]:

# Substituting espaces and dashes in the column names by underscores

data.columns = data.columns.str.replace(" ", "_")

data.columns = data.columns.str.replace("-", "_")

In [0]:

# Checking

data.columns

Out[0]:

Index(['row_id', 'order_id', 'order_date', 'ship_date', 'ship_mode',

       'customer_id', 'customer_name', 'segment', 'country', 'city', 'state',

       'postal_code', 'region', 'product_id', 'category', 'sub_category',

       'product_name', 'sales', 'quantity', 'discount', 'profit'],

      dtype='object')

In [0]:

# Assess the unique values per column (to know whether a variable is categorical or not)

for c in data.columns:

    if len(set(data[c])) < 20:

        print(c,set(data[c]))

ship_mode {'Second Class', 'First Class', 'Same Day', 'Standard Class'}

segment {'Corporate', 'Home Office', 'Consumer'}

country {'United States'}

region {'Central', 'East', 'South', 'West'}

category {'Technology', 'Furniture', 'Office Supplies'}

sub_category {'Fasteners', 'Tables', 'Art', 'Phones', 'Bookcases', 'Machines', 'Copiers', 'Labels', 'Accessories', 'Supplies', 'Appliances', 'Envelopes', 'Paper', 'Furnishings', 'Binders', 'Storage', 'Chairs'}

quantity {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14}

discount {0.0, 0.8, 0.2, 0.3, 0.45, 0.5, 0.7, 0.6, 0.32, 0.1, 0.4, 0.15}

In [0]:

# Splitting data by category

df_technology = data.loc[data['category'] == 'Technology']

df_furniture = data.loc[data['category'] == 'Furniture']

df_office = data.loc[data['category'] == 'Office Supplies']

In [0]:

# Aggregating sales by order date

ts_technology = df_technology.groupby('order_date')['sales'].sum().reset_index()

ts_furniture = df_furniture.groupby('order_date')['sales'].sum().reset_index()

ts_office = df_office.groupby('order_date')['sales'].sum().reset_index()

In [0]:

# Checking dataset

ts_technology

Out[0]:

824 rows × 2 columns

In [0]:

# Setting the date as index

ts_technology = ts_technology.set_index('order_date')

ts_furniture = ts_furniture.set_index('order_date')

ts_office = ts_office.set_index('order_date')

In [0]:

# Visualizing the series

ts_technology

Out[0]:

824 rows × 1 columns

In [0]:

# Plotting sales in the technology category

sales_technology = ts_technology[['sales']]

ax = sales_technology.plot(color = 'b', figsize = (18,6))

plt.xlabel("Data")

plt.ylabel('Sales')

plt.title("Technology Category Sales")

plt.show()

In [0]:

# Plotting sales in the furniture category

sales_furniture = ts_furniture[['sales']]

ax = sales_furniture.plot(color = 'g', figsize = (18,6))

plt.xlabel("Data")

plt.ylabel('Sales')

plt.title("Furniture Category Sales")

plt.show()

In [0]:

# Plotting sales in the office category

sales_office = ts_office[['sales']]

ax = sales_office.plot(color = 'r', figsize = (18,6))

plt.xlabel("Data")

plt.ylabel('Sales')

plt.title("Office Category Sales")

plt.show()

Adjusting the index type to DateTimeIndex (which characterizes a time series), so that it's possible to aggregate monthly and obtain the mean monthly sales.

In [0]:

# Checking index type

type(sales_technology.index)

Out[0]:

pandas.core.indexes.base.Index

In [0]:

# Changing index type

sales_technology.index = pd.to_datetime(sales_technology.index)

sales_furniture.index = pd.to_datetime(sales_furniture.index)

sales_office.index = pd.to_datetime(sales_office.index)

In [0]:

# Checking index type

type(sales_technology.index)

Out[0]:

pandas.core.indexes.datetimes.DatetimeIndex

In [0]:

# Resampling data to a monthly frequency

# Done by setting the month as index and then calculating the mean of daily sales over the month

sales_technology_monthly_mean = sales_technology['sales'].resample('MS').mean()

sales_furniture_monthly_mean = sales_furniture['sales'].resample('MS').mean()

sales_office_monthly_mean = sales_office['sales'].resample('MS').mean()

In [0]:

# Verifying the resulting type

type(sales_technology_monthly_mean)

Out[0]:

pandas.core.series.Series

In [0]:

# Checking the data

sales_technology_monthly_mean

Out[0]:

order_date

2014-01-01     449.041429

2014-02-01     229.787143

2014-03-01    2031.948375

2014-04-01     613.028933

2014-05-01     564.698588

2014-06-01     766.905909

2014-07-01     533.608933

2014-08-01     708.435385

2014-09-01    2035.838133

2014-10-01     596.900900

2014-11-01    1208.056320

2014-12-01    1160.732889

2015-01-01     925.070800

2015-02-01     431.121250

2015-03-01     574.662333

2015-04-01     697.559500

2015-05-01     831.642857

2015-06-01     429.024400

2015-07-01     691.397733

2015-08-01    1108.902286

2015-09-01     950.856400

2015-10-01     594.716111

2015-11-01    1037.982652

2015-12-01    1619.637636

2016-01-01     374.671067

2016-02-01    1225.891400

2016-03-01    1135.150105

2016-04-01     875.911882

2016-05-01    1601.816167

2016-06-01    1023.259500

2016-07-01     829.312500

2016-08-01     483.620100

2016-09-01    1144.170300

2016-10-01    1970.835875

2016-11-01    1085.642360

2016-12-01     970.554870

2017-01-01    1195.218071

2017-02-01     430.501714

2017-03-01    1392.859250

2017-04-01     825.559133

2017-05-01     678.329400

2017-06-01     853.055000

2017-07-01    1054.996636

2017-08-01     978.842333

2017-09-01    1077.704120

2017-10-01    1493.439227

2017-11-01    1996.750920

2017-12-01     955.865652

Freq: MS, Name: sales, dtype: float64

In [0]:

# Plotting the monthly mean daily sales of technology products

sales_technology_monthly_mean.plot(figsize = (18, 6), color = 'blue')

plt.xlabel("Purchase Date")

plt.ylabel('Sales')

plt.title("Technology: Average Daily Sales per Month")

plt.show()

In [0]:

# Plotting the monthly mean daily sales of furniture products

sales_furniture_monthly_mean.plot(figsize = (18, 6), color = 'green')

plt.xlabel("Purchase Date")

plt.ylabel('Sales')

plt.title("Furniture: Average Daily Sales per Month")

plt.show()

In [0]:

# Plotting the monthly mean daily sales of office products

sales_office_monthly_mean.plot(figsize = (20, 6), color = 'red')

plt.xlabel("Purchase Date")

plt.ylabel('Sales')

plt.title("Office: Average Daily Sales per Month")

plt.show()

Decomposing the series to analyze its componentes.

In [0]:

# Decomposing the time series with the monthly mean daily sales of technology products

decomposition = seasonal_decompose(sales_technology_monthly_mean, freq = 12)

rcParams['figure.figsize'] = 18, 12

# Components

trend = decomposition.trend

seasonal = decomposition.seasonal

residual = decomposition.resid

# Plot

plt.subplot(411)

plt.plot(sales_technology_monthly_mean, label = 'Original Series')

plt.legend(loc = 'best')

plt.subplot(412)

plt.plot(trend, label = 'Trend')

plt.legend(loc = 'best')

plt.subplot(413)

plt.plot(seasonal, label = 'Seasonality')

plt.legend(loc = 'best')

plt.subplot(414)

plt.plot(residual, label = 'Residuals')

plt.legend(loc = 'best')

plt.tight_layout()

Stationarity test:

In [0]:

# Function to test stationarity

def stationarity_test(serie):

    # Calculating moving statistics

    rolmean = serie.rolling(window = 12).mean()

    rolstd = serie.rolling(window = 12).std()

    # Plot of moving statistics

    orig = plt.plot(serie, color = 'blue', label = 'Original')

    mean = plt.plot(rolmean, color = 'red', label = 'Moving Average')

    std = plt.plot(rolstd, color = 'black', label = 'Standard Deviation')

    plt.legend(loc = 'best')

    plt.title('Moving Statistics - Mean and Standard Deviation')

    plt.show()

    # Dickey-Fuller test:

    # Print

    print('\nDickey-Fuller Test Result:\n')

    # Test

    dftest = adfuller(serie, autolag = 'AIC')

    # Formatting the output

    df_output = pd.Series(dftest[0:4], index = ['Test Statistic',

                                               'P-value',

                                               'Number of Lags Considered',

                                               'Number of Observations Used'])

    # Loop through each item in the test output

    for key, value in dftest[4].items():

        df_output['Critical Value (%s)'%key] = value

    # Print

    print (df_output)

    # Testing p-value

    print ('\nConclusion:')

    if df_output[1] > 0.05:

        print('\nThe p-value is above 0.05, therefore there are no evidences to reject the null hypothesis.')

        print('This series probably is not stationary.')     

    else:

        print('\nThe p-value is below 0.05, therefore there are evidences to reject the null hypothesis.')

        print('This series probably is stationary.')        

In [0]:

# Verifying if the series is stationary

stationarity_test(sales_technology_monthly_mean)

Dickey-Fuller Test Result:

Test Statistic                -7.187969e+00

P-value                        2.547334e-10

Number of Lags Considered      0.000000e+00

Number of Observations Used    4.700000e+01

Critical Value (1%)           -3.577848e+00

Critical Value (5%)           -2.925338e+00

Critical Value (10%)          -2.600774e+00

dtype: float64

Conclusion:

The p-value is below 0.05, therefore there are evidences to reject the null hypothesis.

This series probably is stationary.

Function to Calculate Accuracy¶

In [0]:

# Function

def performance(y_true, y_pred):

    mse = ((y_pred - y_true) ** 2).mean()

    mape = np.mean(np.abs((y_true - y_pred) / y_true)) * 100

    return( print('The prediction MSE is {}'.format(round(mse, 2))+

                  '\nThe prediction RMSE is {}'.format(round(np.sqrt(mse), 2))+

                  '\nThe prediction MAPE is {}'.format(round(mape, 2))))

ARMA Model¶

In [0]:

# Load data

data = pd.read_csv('https://raw.githubusercontent.com/dsacademybr/Datasets/master/dataset6.csv')

# Setting column names to lowercase

data.columns = map(str.lower, data.columns)

# Substituting espaces and dashes in the column names by underscores

data.columns = data.columns.str.replace(" ", "_")

data.columns = data.columns.str.replace("-", "_")

# Splitting data by category

df_technology = data.loc[data['category'] == 'Technology']

# Aggregating sales by order date

ts_technology = df_technology.groupby('order_date')['sales'].sum().reset_index()

# Setting the date as index

ts_technology = ts_technology.set_index('order_date')

# Retrieving only sales data

sales_technology = ts_technology[['sales']]

# Changing index type

sales_technology.index = pd.to_datetime(sales_technology.index)

# Resampling data to a monthly frequency

# Done by setting the month as index and then calculating the mean of daily sales over the month

sales_technology_monthly_mean = sales_technology['sales'].resample('MS').mean()

# Train-Test Split

X = sales_technology_monthly_mean

train_size = int(len(X) * 0.75)

trainset, testset = X[0:train_size], X[train_size:]

https://www.statsmodels.org/stable/generated/statsmodels.tsa.arima_model.ARMA.html

In [0]:

# Seeking the best values for the 'p' and 'q' parameters

# The goal is the lowest possible AIC

# Initial values

best_aic = np.inf

best_order = None

best_model = None

# Values to be used in seeking the best combination, order = (i, j)

test_order_values = [1, 2, 3, 4, 5]

# Loop

for i in test_order_values:

    for j in test_order_values:

        try:

            tmp_mdl = ARMA(trainset, order = (i, j)).fit(disp = False)

            tmp_aic = tmp_mdl.aic

            if tmp_aic < best_aic:

                best_aic = tmp_aic

                best_order = (i, j)

                best_model = tmp_mdl

        except: continue

# Print

print('\nBest AIC Score: %6.2f | order: %s'%(best_aic, best_order))

print('\nBest Model:', best_model.summary())

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:695: RuntimeWarning:

divide by zero encountered in log

/usr/local/lib/python3.6/dist-packages/statsmodels/tools/numdiff.py:243: RuntimeWarning:

invalid value encountered in subtract

/usr/local/lib/python3.6/dist-packages/statsmodels/tools/numdiff.py:243: RuntimeWarning:

invalid value encountered in multiply

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:668: RuntimeWarning:

invalid value encountered in exp

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:669: RuntimeWarning:

invalid value encountered in exp

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:512: ConvergenceWarning:

Maximum Likelihood optimization failed to converge. Check mle_retvals

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:693: RuntimeWarning:

invalid value encountered in double_scalars

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:668: RuntimeWarning:

overflow encountered in exp

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:668: RuntimeWarning:

invalid value encountered in true_divide

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:669: RuntimeWarning:

overflow encountered in exp

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:669: RuntimeWarning:

invalid value encountered in true_divide

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:695: RuntimeWarning:

divide by zero encountered in true_divide

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:512: ConvergenceWarning:

Maximum Likelihood optimization failed to converge. Check mle_retvals

Best AIC Score: 552.92 | order: (3, 1)

Best Model:                               ARMA Model Results                             

==============================================================================

Dep. Variable:                  sales   No. Observations:                   36

Model:                     ARMA(3, 1)   Log Likelihood                -270.459

Method:                       css-mle   S.D. of innovations            442.123

Date:                Sat, 09 May 2020   AIC                            552.918

Time:                        17:27:13   BIC                            562.420

Sample:                    01-01-2014   HQIC                           556.235

                         - 12-01-2016                                        

===============================================================================

                  coef    std err          z      P>|z|      [0.025      0.975]

-------------------------------------------------------------------------------

const         939.5686     63.228     14.860      0.000     815.644    1063.493

ar.L1.sales    -0.3149      0.770     -0.409      0.685      -1.824       1.194

ar.L2.sales    -0.2023      0.207     -0.979      0.335      -0.607       0.203

ar.L3.sales     0.1073      0.260      0.413      0.682      -0.402       0.616

ma.L1.sales     0.1991      0.759      0.262      0.795      -1.289       1.687

                                    Roots                                   

=============================================================================

                  Real          Imaginary           Modulus         Frequency

-----------------------------------------------------------------------------

AR.1           -0.8027           -1.4232j            1.6340           -0.3317

AR.2           -0.8027           +1.4232j            1.6340            0.3317

AR.3            3.4918           -0.0000j            3.4918           -0.0000

MA.1           -5.0226           +0.0000j            5.0226            0.5000

-----------------------------------------------------------------------------

The warning messages above indicate taht for some order combinations (i, j) the model could not converge. The combination with the best performance was: order = (3, 1)

In [0]:

# Create model

arma_model = ARMA(trainset, order = (3, 1))

In [0]:

# Train model

arma_model_fit = arma_model.fit(method = 'mle', disp = False)

In [0]:

# Model summary

arma_model_fit.summary()

Out[0]:

Evaluating based on the AIC (Akaike Information Criterion) metric. The lower the value the better the model performance.

In [0]:

# Model predictions

arma_predict = arma_model_fit.predict(start = pd.to_datetime('2017-01-01'),

                                       end = pd.to_datetime('2017-12-01'),

                                       dynamic = False)

In [0]:

# Plotting predictions

ax = sales_technology_monthly_mean.plot(label = 'Observed Values', color = '#3574BF')

rcParams['figure.figsize'] = 18, 8

arma_predict.plot(ax = ax, label = 'ARMA (3,1) Predictions', alpha = 0.7, color = 'red')

plt.title('Sales Predictions with ARMA Model')

plt.xlabel('Data')

plt.ylabel('Technology Category Sales')

plt.legend()

plt.show()

In [0]:

# Evaluating performance with test data

arma_results = performance(testset, arma_predict)

arma_results

The prediction MSE is 168192.9

The prediction RMSE is 410.11

The prediction MAPE is 28.4

Creating a function to plot various diagnostic metrics which can be used to assess the model

In [0]:

# Function

def tsplot(y, lags = None, figsize = (12, 8), style = 'bmh'):

    # If the series isn't from type pd.Series, convert it

    if not isinstance(y, pd.Series):

        y = pd.Series(y)

    # Creating plots

    with plt.style.context(style):   

        fig = plt.figure(figsize = figsize)

        layout = (3, 2)

        ts_ax = plt.subplot2grid(layout, (0, 0), colspan = 2)

        acf_ax = plt.subplot2grid(layout, (1, 0))

        pacf_ax = plt.subplot2grid(layout, (1, 1))

        qq_ax = plt.subplot2grid(layout, (2, 0))

        pp_ax = plt.subplot2grid(layout, (2, 1))

        y.plot(ax = ts_ax)

        ts_ax.set_title('Plots for Time Series Analysis')

        smt.graphics.plot_acf(y, lags = lags, ax = acf_ax, alpha = 0.05)

        smt.graphics.plot_pacf(y, lags = lags, ax = pacf_ax, alpha = 0.05)

        sm.qqplot(y, line = 's', ax = qq_ax)

        qq_ax.set_title('QQ Plot')       

        scs.probplot(y, sparams = (y.mean(), y.std()), plot = pp_ax)

        plt.tight_layout()

    return

In [0]:

# Running the functions with train data

tsplot(trainset, lags = 30)

/usr/local/lib/python3.6/dist-packages/statsmodels/regression/linear_model.py:1358: RuntimeWarning:

invalid value encountered in sqrt

The model's rediduals must be normally distributed. Checking with the jarque_bera test.

In [0]:

# Test

score, pvalue, _, _ = jarque_bera(arma_model_fit.resid)

# Result

if pvalue < 0.05:

    print ('\nThe residuals might not be normally distributed.')

else:

    print ('\nThe residuals appear to be normally distributed.') 

The residuals appear to be normally distributed.

In [0]:

# Ljung-Box test

ljungbox_test = sms.diagnostic.acorr_ljungbox(arma_model_fit.resid, lags = [30], boxpierce = False)

print('P-value =', ljungbox_test[1])

P-value = [0.9105266]

In [0]:

# Load data

data = pd.read_csv('https://raw.githubusercontent.com/dsacademybr/Datasets/master/dataset6.csv')

# Setting column names to lowercase

data.columns = map(str.lower, data.columns)

# Substituting espaces and dashes in the column names by underscores

data.columns = data.columns.str.replace(" ", "_")

data.columns = data.columns.str.replace("-", "_")

# Splitting data by category

df_technology = data.loc[data['category'] == 'Technology']

# Aggregating sales by order date

ts_technology = df_technology.groupby('order_date')['sales'].sum().reset_index()

# Setting the date as index

ts_technology = ts_technology.set_index('order_date')

# Retrieving only sales data

sales_technology = ts_technology[['sales']]

# Changing index type

sales_technology.index = pd.to_datetime(sales_technology.index)

# Resampling data to a monthly frequency

# Done by setting the month as index and then calculating the mean of daily sales over the month

sales_technology_monthly_mean = sales_technology['sales'].resample('MS').mean()

# Train-Test Split

X = sales_technology_monthly_mean

train_size = int(len(X) * 0.75)

trainset, testset = X[0:train_size], X[train_size:]

In [0]:

# Function to evaluate an ARIMA model

def evaluate_arima_model(X, arima_order):

    # Prepare data

    train_size = int(len(X) * 0.75)

    train, test = X[0:train_size], X[train_size:]

    history = [x for x in train]

    # predict

    predictions = list()

    # Loop

    for t in range(len(test)):

        model = ARIMA(history, order = arima_order)

        model_fit = model.fit(method = 'mle', disp = 0)

        yhat = model_fit.forecast()[0]

        predictions.append(yhat)

        history.append(test[t])

    # Calculate model error

    error = mean_squared_error(test, predictions)

    return error

In [0]:

# Finding the optimal combination of 'p', 'd' and 'q' values for the ARIMA model

def find_best_arima(dataset, p_values, d_values, q_values):

    # Adjusting the data type

    dataset = dataset.astype('float32')

    # Defining control variables

    best_score, best_cfg = float("inf"), None

    # Loop through 'p d q' values

    for p in p_values:

        for d in d_values:

            for q in q_values:

                order = (p,d,q)

                try:

                    mse = evaluate_arima_model(dataset, order)

                    if mse < best_score:

                        best_score, best_cfg = mse, order

                    print('ARIMA%s MSE = %.3f' % (order, mse))

                except:

                    continue

    print('\nBest Model ARIMA%s MSE = %.3f' % (best_cfg, best_score))

In [0]:

# Set values to loop through

p_values = [0, 2, 4, 6]

d_values = [0, 2, 4, 6]

q_values = [0, 2, 4, 6] 

In [0]:

# Applying function to seek the best values for 'p d q'

find_best_arima(sales_technology_monthly_mean.values, p_values, d_values, q_values)

ARIMA(0, 0, 0) MSE = 174155.275

ARIMA(0, 0, 2) MSE = 191647.549

ARIMA(0, 0, 4) MSE = 194995.823

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

ARIMA(0, 0, 6) MSE = 173741.615

ARIMA(0, 2, 0) MSE = 793878.683

ARIMA(2, 0, 0) MSE = 194636.248

ARIMA(2, 2, 0) MSE = 299872.253

ARIMA(4, 0, 0) MSE = 188135.157

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

ARIMA(4, 0, 2) MSE = 201385.415

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:668: RuntimeWarning:

overflow encountered in exp

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:668: RuntimeWarning:

invalid value encountered in true_divide

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:669: RuntimeWarning:

overflow encountered in exp

/usr/local/lib/python3.6/dist-packages/statsmodels/tsa/tsatools.py:669: RuntimeWarning:

invalid value encountered in true_divide

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

ARIMA(4, 0, 4) MSE = 205003.913

ARIMA(4, 2, 0) MSE = 299179.994

ARIMA(6, 0, 0) MSE = 204155.758

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

ARIMA(6, 0, 2) MSE = 183004.507

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

ARIMA(6, 0, 4) MSE = 140567.260

Best Model ARIMA(6, 0, 4) MSE = 140567.260

/usr/local/lib/python3.6/dist-packages/statsmodels/base/model.py:492: HessianInversionWarning:

Inverting hessian failed, no bse or cov_params available

In [0]:

# Create model

arima_model = ARIMA(trainset, order = (6,0,4))

In [0]:

# Train model

arima_model_fit = arima_model.fit(method = 'mle', disp = False)

In [0]:

# Model summary

arima_model_fit.summary()

Out[0]:

Evaluating based on the AIC (Akaike Information Criterion) metric. The lower the value the better the model performance.

In [0]:

# Model predictions

arima_predict = arima_model_fit.predict(start = pd.to_datetime('2017-01-01'),

                                         end = pd.to_datetime('2017-12-01'),

                                         dynamic = False)

In [0]:

# Plotting predictions

ax = sales_technology_monthly_mean.plot(label = 'Observed Values', color = '#3574BF')

rcParams['figure.figsize'] = 14, 7

arima_predict.plot(ax = ax, label = 'ARIMA (6,0,4) Predictions', alpha = 0.8, color = 'red')

plt.title('Sales Predictions with ARIMA Model')

plt.xlabel('Data')

plt.ylabel('Technology Category Sales')

plt.legend()

plt.show()

In [0]:

# Calculating performance

arima_results = performance(testset, arima_predict)

arima_results

The prediction MSE is 120137.97

The prediction RMSE is 346.61

The prediction MAPE is 20.1

The model's rediduals must be normally distributed. Checking with the jarque_bera test..

In [0]:

# Test

score, pvalue, _, _ = jarque_bera(arima_model_fit.resid)

# Result

if pvalue < 0.05:

    print ('\nThe residuals might not be normally distributed.')

else:

    print ('\nThe residuals appear to be normally distributed.')

The residuals appear to be normally distributed.

In [0]:

# Ljung-Box test

ljungbox_test = sms.diagnostic.acorr_ljungbox(arima_model_fit.resid, lags = [30], boxpierce = False)

print('P-value =', ljungbox_test[1])

P-value = [0.92902243]

ARMA (3,1):

• The prediction MSE is 168192.9
• The prediction RMSE is 410.11
• The prediction MAPE is 28.4

ARIMA (6,0,4):

• The prediction MSE is 120137.97
• The prediction RMSE is 346.61
• The prediction MAPE is 20.1

In [0]:

# Load data

data = pd.read_csv('https://raw.githubusercontent.com/dsacademybr/Datasets/master/dataset6.csv')

# Setting column names to lowercase

data.columns = map(str.lower, data.columns)

# Substituting espaces and dashes in the column names by underscores

data.columns = data.columns.str.replace(" ", "_")

data.columns = data.columns.str.replace("-", "_")

# Splitting data by category

df_technology = data.loc[data['category'] == 'Technology']

# Aggregating sales by order date

ts_technology = df_technology.groupby('order_date')['sales'].sum().reset_index()

# Setting the date as index

ts_technology = ts_technology.set_index('order_date')

# Retrieving only sales data

sales_technology = ts_technology[['sales']]

# Changing index type

sales_technology.index = pd.to_datetime(sales_technology.index)

# Resampling data to a monthly frequency

# Done by setting the month as index and then calculating the mean of daily sales over the month

sales_technology_monthly_mean = sales_technology['sales'].resample('MS').mean()

# Train-Test Split

X = sales_technology_monthly_mean

train_size = int(len(X) * 0.75)

trainset, testset = X[0:train_size], X[train_size:]

In [0]:

# Searching for the ideal model order

# Using pm.auto_arima method to apply a Grid Search and retrieve the best model

sarima_hiperparameters = pm.auto_arima(trainset,

                          seasonal = True,

                          m = 12,

                          d = 0,

                          D = 1,

                          max_p = 2,

                          max_q = 2,

                          trace = True,

                          error_action = 'ignore',

                          suppress_warnings = True)

Performing stepwise search to minimize aic

Fit ARIMA(2,0,2)x(1,1,1,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(0,0,0)x(0,1,0,12) [intercept=True]; AIC=380.244, BIC=382.600, Time=0.024 seconds

Fit ARIMA(1,0,0)x(1,1,0,12) [intercept=True]; AIC=378.184, BIC=382.897, Time=0.495 seconds

Fit ARIMA(0,0,1)x(0,1,1,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(0,0,0)x(0,1,0,12) [intercept=False]; AIC=378.609, BIC=379.787, Time=0.015 seconds

Fit ARIMA(1,0,0)x(0,1,0,12) [intercept=True]; AIC=382.218, BIC=385.753, Time=0.068 seconds

Fit ARIMA(1,0,0)x(2,1,0,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(1,0,0)x(1,1,1,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(1,0,0)x(0,1,1,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(1,0,0)x(2,1,1,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(0,0,0)x(1,1,0,12) [intercept=True]; AIC=377.369, BIC=380.903, Time=0.331 seconds

Fit ARIMA(0,0,0)x(2,1,0,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(0,0,0)x(1,1,1,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(0,0,0)x(0,1,1,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(0,0,0)x(2,1,1,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(0,0,1)x(1,1,0,12) [intercept=True]; AIC=376.797, BIC=381.509, Time=0.255 seconds

Fit ARIMA(0,0,1)x(0,1,0,12) [intercept=True]; AIC=382.184, BIC=385.718, Time=0.145 seconds

Fit ARIMA(0,0,1)x(2,1,0,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(0,0,1)x(1,1,1,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(0,0,1)x(2,1,1,12) [intercept=True]; AIC=nan, BIC=nan, Time=nan seconds

Fit ARIMA(1,0,1)x(1,1,0,12) [intercept=True]; AIC=378.567, BIC=384.458, Time=0.472 seconds

Fit ARIMA(0,0,2)x(1,1,0,12) [intercept=True]; AIC=378.319, BIC=384.209, Time=0.330 seconds

Fit ARIMA(1,0,2)x(1,1,0,12) [intercept=True]; AIC=379.996, BIC=387.065, Time=0.518 seconds

Total fit time: 2.702 seconds

Analyzing the results from the Stepwise Grid Search the best model is:

Fit ARIMA(0,0,1)x(1,1,0,12) [intercept=True]; AIC=376.797, BIC=381.509, Time=0.225 seconds

In [0]:

# Print summary of the model with best parameters

print(sarima_hiperparameters.summary())

                                 Statespace Model Results                                

==========================================================================================

Dep. Variable:                                  y   No. Observations:                   36

Model:             SARIMAX(0, 0, 1)x(1, 1, 0, 12)   Log Likelihood                -184.398

Date:                            Sat, 09 May 2020   AIC                            376.797

Time:                                    17:29:14   BIC                            381.509

Sample:                                         0   HQIC                           378.047

                                             - 36                                        

Covariance Type:                              opg                                        

==============================================================================

                 coef    std err          z      P>|z|      [0.025      0.975]

------------------------------------------------------------------------------

intercept    163.5162     74.251      2.202      0.028      17.986     309.046

ma.L1         -0.4516      0.232     -1.948      0.051      -0.906       0.003

ar.S.L12      -0.7046      0.146     -4.812      0.000      -0.992      -0.418

sigma2      2.449e+05   1.09e+05      2.253      0.024    3.19e+04    4.58e+05

===================================================================================

Ljung-Box (Q):                         nan   Jarque-Bera (JB):                 0.65

Prob(Q):                               nan   Prob(JB):                         0.72

Heteroskedasticity (H):               2.44   Skew:                             0.13

Prob(H) (two-sided):                  0.23   Kurtosis:                         2.24

===================================================================================

Warnings:

[1] Covariance matrix calculated using the outer product of gradients (complex-step).

https://www.statsmodels.org/dev/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html

In [0]:

# Create SARIMA model

sarima_model_v1 = sm.tsa.statespace.SARIMAX(trainset,

                                             order = (0, 0, 1),

                                             seasonal_order = (1, 1, 0, 12),

                                             enforce_stationarity = False,

                                             enforce_invertibility = False)

In [0]:

# Training (fit) model

sarima_model_v1_fit = sarima_model_v1.fit()

In [0]:

# Model summary

print(sarima_model_v1_fit.summary())

                                 Statespace Model Results                                

==========================================================================================

Dep. Variable:                              sales   No. Observations:                   36

Model:             SARIMAX(0, 0, 1)x(1, 1, 0, 12)   Log Likelihood                 -93.817

Date:                            Sat, 09 May 2020   AIC                            193.633

Time:                                    17:29:14   BIC                            195.088

Sample:                                01-01-2014   HQIC                           193.095

                                     - 12-01-2016                                        

Covariance Type:                              opg                                        

==============================================================================

                 coef    std err          z      P>|z|      [0.025      0.975]

------------------------------------------------------------------------------

ma.L1         -0.0756      0.524     -0.144      0.885      -1.103       0.952

ar.S.L12      -0.4330      0.602     -0.719      0.472      -1.613       0.747

sigma2      3.385e+05   1.55e+05      2.183      0.029    3.46e+04    6.42e+05

===================================================================================

Ljung-Box (Q):                         nan   Jarque-Bera (JB):                 0.93

Prob(Q):                               nan   Prob(JB):                         0.63

Heteroskedasticity (H):               2.43   Skew:                             0.55

Prob(H) (two-sided):                  0.41   Kurtosis:                         2.20

===================================================================================

Warnings:

[1] Covariance matrix calculated using the outer product of gradients (complex-step).

In [0]:

# Model diagnostic

sarima_model_v1_fit.plot_diagnostics(figsize = (16, 8))

plt.show()

In [0]:

# Predictions

sarima_predict_1 = sarima_model_v1_fit.get_prediction(start = pd.to_datetime('2017-01-01'),

                                                       end = pd.to_datetime('2017-12-01'),

                                                       dynamic = False)

In [0]:

# Confidence interval

sarima_predict_conf_1 = sarima_predict_1.conf_int()

sarima_predict_conf_1

Out[0]:

In [0]:

# Plot observed values

ax = sales_technology_monthly_mean.plot(label = 'Observed Values', color = '#2574BF')

# Plot predicted values

sarima_predict_1.predicted_mean.plot(ax = ax,

                                     label = 'SARIMA (0, 0, 1)x(0, 1, 1, 12) Predictions',

                                     alpha = 0.7,

                                     color = 'red')

# Plot confidence interval

ax.fill_between(sarima_predict_conf_1.index,

                # lower sales

                sarima_predict_conf_1.iloc[:, 0],

                # upper sales

                sarima_predict_conf_1.iloc[:, 1], color = 'k', alpha = 0.1)

# Titles and legends

plt.title('Sales Predictions with SARIMA Model')

plt.xlabel('Data')

plt.ylabel('Technology Category Sales')

plt.legend()

plt.show()

In [0]:

# Calculating performance

sarima_results_1 = performance(testset, sarima_predict_1.predicted_mean)

sarima_results_1

The prediction MSE is 184433.09

The prediction RMSE is 429.46

The prediction MAPE is 35.35

In [0]:

# Ljung-Box test

ljungbox_test = sms.diagnostic.acorr_ljungbox(sarima_model_v1_fit.resid, lags = [30], boxpierce = False)

print('P-value =', ljungbox_test[1])

P-value = [0.6456201]

ARMA (3,1):

• The prediction MSE is 168192.9
• The prediction RMSE is 410.11
• The prediction MAPE is 28.4

ARIMA (6,0,4):

• The prediction MSE is 120137.97
• The prediction RMSE is 346.61
• The prediction MAPE is 20.1

SARIMA (0,0,1)x(1,1,0,12):

• The prediction MSE is 184433.09
• The prediction RMSE is 429.46
• The prediction MAPE is 35.35

Grid Search Method 2 - manual¶

In [0]:

# Load data

data = pd.read_csv('https://raw.githubusercontent.com/dsacademybr/Datasets/master/dataset6.csv')

# Setting column names to lowercase

data.columns = map(str.lower, data.columns)

# Substituting espaces and dashes in the column names by underscores

data.columns = data.columns.str.replace(" ", "_")

data.columns = data.columns.str.replace("-", "_")

# Splitting data by category

df_technology = data.loc[data['category'] == 'Technology']

# Aggregating sales by order date

ts_technology = df_technology.groupby('order_date')['sales'].sum().reset_index()

# Setting the date as index

ts_technology = ts_technology.set_index('order_date')

# Retrieving only sales data

sales_technology = ts_technology[['sales']]

# Changing index type

sales_technology.index = pd.to_datetime(sales_technology.index)

# Resampling data to a monthly frequency

# Done by setting the month as index and then calculating the mean of daily sales over the month

sales_technology_monthly_mean = sales_technology['sales'].resample('MS').mean()

# Train-Test Split

X = sales_technology_monthly_mean

train_size = int(len(X) * 0.75)

trainset, testset = X[0:train_size], X[train_size:]

Testing multiple hiperparameter combinations for the Sarima model. Since going through all possible hiperparameter combinations for all variables in a sarima model can take long, I'm limiting the value of each hiperparameter to either 0 or 1, since that range fits most datasets.

In [0]:

# Setting 'p d q' to be either 0 or 1 (range(0, 2))

# Since there are many parameters (p,d,q)(P,D,Q) anything more than this start taking a long time to train

p = d = q = range(0, 2)

In [0]:

# List with all combinations for 'p d q'

pdq = list(itertools.product(p, d, q))

pdq

Out[0]:

[(0, 0, 0),

 (0, 0, 1),

 (0, 1, 0),

 (0, 1, 1),

 (1, 0, 0),

 (1, 0, 1),

 (1, 1, 0),

 (1, 1, 1)]

In [0]:

# List with all combinations for the seasonal hiperparameters 'P D Q'

# Using List Comprehension

seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]

seasonal_pdq

Out[0]:

[(0, 0, 0, 12),

 (0, 0, 1, 12),

 (0, 1, 0, 12),

 (0, 1, 1, 12),

 (1, 0, 0, 12),

 (1, 0, 1, 12),

 (1, 1, 0, 12),

 (1, 1, 1, 12)]

In [0]:

# Grid Search

#warnings.filterwarnings("ignore")

# Lowest possible value for the AIC statistics (the model optimization goal)

lowest_aic = sys.maxsize

lowest = ''

# Loop

for param in pdq:

    for param_seasonal in seasonal_pdq:

        try:

            # Create model with the hiperparameter combination

            mod = sm.tsa.statespace.SARIMAX(trainset,

                                            order = param,

                                            seasonal_order = param_seasonal,

                                            enforce_stationarity = False,

                                            enforce_invertibility = False)

            # Train model

            results = mod.fit()

            # Print

            print('SARIMA{}x{} - AIC:{}'.format(param, param_seasonal, results.aic))

            # Check for improvement (reduction) in AIC

            if lowest_aic >  results.aic:

                lowest = 'SARIMA{}x{} - AIC:{}'.format(param, param_seasonal, results.aic)

                lowest_aic = results.aic

        except:

            continue

print ("\nModel with Lowest AIC Score: " + lowest)

SARIMA(0, 0, 0)x(0, 0, 0, 12) - AIC:588.2102390405763

SARIMA(0, 0, 0)x(0, 1, 0, 12) - AIC:363.3206300182339

SARIMA(0, 0, 0)x(1, 0, 0, 12) - AIC:379.0151254818642

SARIMA(0, 0, 0)x(1, 1, 0, 12) - AIC:191.64720213115754

SARIMA(0, 0, 1)x(0, 0, 0, 12) - AIC:556.7477551380659

SARIMA(0, 0, 1)x(0, 1, 0, 12) - AIC:350.40993759684073

SARIMA(0, 0, 1)x(1, 0, 0, 12) - AIC:379.8551441468676

SARIMA(0, 0, 1)x(1, 1, 0, 12) - AIC:193.63316082077588

SARIMA(0, 1, 0)x(0, 0, 0, 12) - AIC:543.2769140078082

SARIMA(0, 1, 0)x(0, 1, 0, 12) - AIC:364.150186670076

SARIMA(0, 1, 0)x(1, 0, 0, 12) - AIC:358.5565559049278

SARIMA(0, 1, 0)x(1, 1, 0, 12) - AIC:182.84893834777014

SARIMA(0, 1, 1)x(0, 0, 0, 12) - AIC:497.7775888290871

SARIMA(0, 1, 1)x(0, 1, 0, 12) - AIC:331.55111404956654

SARIMA(0, 1, 1)x(1, 0, 0, 12) - AIC:348.1154222267687

SARIMA(0, 1, 1)x(1, 1, 0, 12) - AIC:178.52195610425457

SARIMA(1, 0, 0)x(0, 0, 0, 12) - AIC:556.6431138617959

SARIMA(1, 0, 0)x(0, 1, 0, 12) - AIC:365.3000316836825

SARIMA(1, 0, 0)x(1, 0, 0, 12) - AIC:358.9225125802152

SARIMA(1, 0, 0)x(1, 1, 0, 12) - AIC:178.59201216997738

SARIMA(1, 0, 1)x(0, 0, 0, 12) - AIC:522.1371166193512

SARIMA(1, 0, 1)x(0, 1, 0, 12) - AIC:352.33098411264257

SARIMA(1, 0, 1)x(1, 0, 0, 12) - AIC:349.64405196314834

SARIMA(1, 0, 1)x(1, 1, 0, 12) - AIC:180.30466198764935

SARIMA(1, 1, 0)x(0, 0, 0, 12) - AIC:535.9247781811981

SARIMA(1, 1, 0)x(0, 1, 0, 12) - AIC:362.75219388274553

SARIMA(1, 1, 0)x(1, 0, 0, 12) - AIC:341.2622085096664

SARIMA(1, 1, 0)x(1, 1, 0, 12) - AIC:166.69418783304513

SARIMA(1, 1, 1)x(0, 0, 0, 12) - AIC:499.3069548304077

SARIMA(1, 1, 1)x(0, 1, 0, 12) - AIC:333.53654481308536

SARIMA(1, 1, 1)x(1, 0, 0, 12) - AIC:334.45693269725643

SARIMA(1, 1, 1)x(1, 1, 0, 12) - AIC:163.8458481947634

Model with Lowest AIC Score: SARIMA(1, 1, 1)x(1, 1, 0, 12) - AIC:163.8458481947634

In [0]:

# Create the model using the best hiperparameters

sarima_model_v2 = sm.tsa.statespace.SARIMAX(trainset,

                                             order = (1, 1, 1),

                                             seasonal_order = (1, 1, 0, 12),

                                             enforce_stationarity = False,

                                             enforce_invertibility = False)

In [0]:

# Train (fit) model

sarima_model_v2_fit = sarima_model_v2.fit()

In [0]:

# Model summary

print(sarima_model_v2_fit.summary())

                                 Statespace Model Results                                

==========================================================================================

Dep. Variable:                              sales   No. Observations:                   36

Model:             SARIMAX(1, 1, 1)x(1, 1, 0, 12)   Log Likelihood                 -77.923

Date:                            Sat, 09 May 2020   AIC                            163.846

Time:                                    17:29:17   BIC                            165.056

Sample:                                01-01-2014   HQIC                           162.518

                                     - 12-01-2016                                        

Covariance Type:                              opg                                        

==============================================================================

                 coef    std err          z      P>|z|      [0.025      0.975]

------------------------------------------------------------------------------

ar.L1          0.0311      0.614      0.051      0.960      -1.172       1.234

ma.L1         -1.0000      0.662     -1.511      0.131      -2.297       0.297

ar.S.L12      -0.2469      0.749     -0.330      0.742      -1.714       1.220

sigma2      3.127e+05   2.12e-06   1.48e+11      0.000    3.13e+05    3.13e+05