<jupyter_start><jupyter_text>Forecasting the monthly sales of French champagne TO-DO* Give a detailed explanation for each step* Fix x axis on line plot* Replace TimeGrouper for Grouper(freq=..) (currently having problems with Grouper)* Additional data exploration* Train a time window based Multilayer Perceptron and compare to ARIMABefore we use ARIMA to model our dataset, let's start by making a naive forecast using walk-forward validation since we want to have a baseline model:<jupyter_code>import scipy
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import sklearn as skl
from sklearn.metrics import mean_squared_error
import statsmodels as stmd
from math import sqrt

def parser(dates):
    return pd.datetime.strptime(dates, '%Y-%m')

series = pd.read_csv('monthly-champagne-sales.csv', header=0, parse_dates=[0], index_col=0, date_parser=parser)
# print(series)
dataset, validation = series[0:-12], series[-12:] #test set is the final year
dataset.to_csv('dataset.csv')
validation.to_csv('validation.csv')

# load data
series = pd.read_csv('dataset.csv', index_col=0, header=0, squeeze=True)
print(series)

# preparing data
X = series.values
# print(X)
X = X.astype('float32')
train_size = int(len(X) * 0.5)
train, test = X[0:train_size], X[train_size:]

# walk-forward validation
history = [x for x in train]
predictions = list()
for i in range(len(test)):
    # predict
    yhat = history[-1]
    predictions.append(yhat)
    # observation
    obs = test[i]
    history.append(obs)
#     print('>Predicted=%.3f, Expected=%.3f' % (yhat, obs))

# report performance
rmse = sqrt(mean_squared_error(test, predictions))
print('RMSE: %.3f' % rmse)<jupyter_output>Month
1964-01-01     2815
1964-02-01     2672
1964-03-01     2755
1964-04-01     2721
1964-05-01     2946
1964-06-01     3036
1964-07-01     2282
1964-08-01     2212
1964-09-01     2922
1964-10-01     4301
1964-11-01     5764
1964-12-01     7312
1965-01-01     2541
1965-02-01     2475
1965-03-01     3031
1965-04-01     3266
1965-05-01     3776
1965-06-01     3230
1965-07-01     3028
1965-08-01     1759
1965-09-01     3595
1965-10-01     4474
1965-11-01     6838
1965-12-01     8357
1966-01-01     3113
1966-02-01     3006
1966-03-01     4047
1966-04-01     3523
1966-05-01     3937
1966-06-01     3986
              ...  
1969-04-01     4276
1969-05-01     4968
1969-06-01     4677
1969-07-01     3523
1969-08-01     1821
1969-09-01     5222
1969-10-01     6872
1969-11-01    10803
1969-12-01    13916
1970-01-01     2639
1970-02-01     2899
1970-03-01     3370
1970-04-01     3740
1970-05-01     2927
1970-06-01     3986
1970-07-01     4217
1970-08-01     1738
1970-09-01     5221
1970-10-01    [...]<jupyter_text>Data Analysis Summary statistics<jupyter_code>print(series.describe())<jupyter_output>count       93.000000
mean      4641.118280
std       2486.403841
min       1573.000000
25%       3036.000000
50%       4016.000000
75%       5048.000000
max      13916.000000
Name: Perrin Freres monthly champagne sales millions ?64-?72, dtype: float64<jupyter_text>Line plot of time seriesLet's look at the monthly sales of Perrin Freres champagne (excluding the data from the validation set):<jupyter_code>series = pd.read_csv('dataset.csv', header=0, squeeze=True)
series.plot()
plt.show()<jupyter_output><empty_output><jupyter_text>It seems there is a seasonality in this time series. *need to fix x-axis  Box and Whisker Plots<jupyter_code>series = pd.read_csv('dataset.csv', header=0, parse_dates=[0], index_col=0, squeeze=True)
# print(series)
groups = series['1964':'1970'].groupby(pd.TimeGrouper('A'))
years = pd.DataFrame()
for name, group in groups:
    years[name.year] = group.values
years.boxplot()
plt.show()<jupyter_output>/Users/jacquesthibodeau/anaconda3/lib/python3.6/site-packages/ipykernel/__main__.py:3: FutureWarning: pd.TimeGrouper is deprecated and will be removed; Please use pd.Grouper(freq=...)
  app.launch_new_instance()<jupyter_text>ARIMAIn order to use ARIMA, we need to make sure we are working with a stationary time series. This is not the case, so we need to difference the series. We'll be differencing the series by substracting the observation from the same time in the previous cycle since it seems like the seasonality is year-to-year. We will lose the first year because it has no prior year to be differenced with, but it is a small price to pay.<jupyter_code>from statsmodels.tsa.stattools import adfuller

# create a differenced series
def difference(dataset, interval=1):
    diff = list()
    for i in range(interval, len(dataset)):
        value = dataset[i] - dataset[i - interval]
        diff.append(value)
    return pd.Series(diff)

series = pd.read_csv('dataset.csv', header=0, parse_dates=[0], index_col=0, squeeze=True)
X = series.values
X = X.astype('float32')

# difference data
months_in_year = 12
stationary = difference(X, months_in_year)
stationary.index = series.index[months_in_year:]
# print(stationary)

# check if stationary
# Dickey-Fuller test for statistical significance
result = adfuller(stationary)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
    print('\t%s: %.3f' % (key, value))
# save
stationary.to_csv('stationary.csv')
# plot
stationary.plot()
plt.show()<jupyter_output>ADF Statistic: -7.134898
p-value: 0.000000
Critical Values:
	1%: -3.515
	5%: -2.898
	10%: -2.586<jupyter_text>Our time series is now stationary. Before we use ARIMA, we need to select the lag values for the Autoregression (p) and Moving Average (q) parameters. With the help of the Autocorrelation Function and the Partial Autocorrelation Function plots, we can get a good starting point.<jupyter_code>from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf

series = pd.read_csv('stationary.csv', header=0, parse_dates=[0], index_col=0)
plot_acf(series, ax=plt.gca())
plt.show()
plot_pacf(series, ax=plt.gca())<jupyter_output><empty_output><jupyter_text>Both figures show significant lag for 1 month. PACF also has some lag at 12 and 13 that might be significant.  ARIMA ModelWe'll be using a grid search to find the best combination of the ARIMA hyperparameters.<jupyter_code>import warnings
from statsmodels.tsa.arima_model import ARIMA
import numpy as np

# create a differenced series
def difference(dataset, interval=1):
    diff = list()
    for i in range(interval, len(dataset)):
        value = dataset[i] - dataset[i - interval]
        diff.append(value)
    return np.array(diff)

# invert differenced value
def inverse_difference(history, yhat, interval=1):
    return yhat + history[-interval]

def arima_model(X, arima_order):
    X = X.astype('float32')
    train_size = int(len(X) * 0.5)
    train, test = X[0:train_size], X[train_size:]

    # walk-forward validation
    history = [x for x in train]
    predictions = list()
    for i in range(len(test)):
        months_in_year = 12
        diff = difference(history, months_in_year)
        # predict
        model = ARIMA(diff, order=arima_order)
        model_fit = model.fit(trend='nc', disp=0)
        yhat = model_fit.forecast()[0]
        yhat = inverse_difference(history, yhat, months_in_year)
        predictions.append(yhat)
        # observation
        obs = test[i]
        history.append(obs)
    #     print('>Predicted=%.3f, Expected=%.3f' % (yhat, obs))

    # report performance
    rmse = sqrt(mean_squared_error(test, predictions))
    return rmse

# evaluating different combinations of the p, d and q values for 
# an ARIMA model's order

def evaluate_arima_models(dataset, p_values, d_values, q_values):
    dataset = dataset.astype('float32')
    best_score, best_cfg = float(1000000), None
    for p in p_values:
        for d in d_values:
            for q in q_values:
                order = (p,d,q)
                try:
                    rmse = arima_model(dataset, order)
                    if rmse < best_score:
                        best_score, best_cfg = rmse, order
                    print('ARIMA%s RMSE=%.3f' % (order,rmse))
                except:
                    continue
    print('Best ARIMA%s RMSE=%.3f' % (best_cfg, best_score))
    
# load data
series = pd.read_csv('dataset.csv', index_col=0)
# print(series)
# evaluate parameters
p_values = range(0, 7)
d_values = range(0, 3)
q_values = range(0, 7)
warnings.filterwarnings("ignore")
evaluate_arima_models(series.values, p_values, d_values, q_values)<jupyter_output>ARIMA(0, 0, 1) RMSE=939.464
ARIMA(0, 0, 2) RMSE=962.289
ARIMA(0, 0, 3) RMSE=944.237
ARIMA(0, 0, 4) RMSE=958.784
ARIMA(0, 0, 5) RMSE=986.832
ARIMA(0, 0, 6) RMSE=1087.037
ARIMA(0, 1, 1) RMSE=958.300
ARIMA(0, 2, 1) RMSE=1146.987
ARIMA(1, 0, 0) RMSE=944.449
ARIMA(1, 1, 0) RMSE=1070.204
ARIMA(1, 1, 1) RMSE=956.945
ARIMA(1, 2, 0) RMSE=1555.858
ARIMA(2, 0, 0) RMSE=955.134
ARIMA(2, 1, 0) RMSE=1030.467
ARIMA(2, 1, 1) RMSE=980.956
ARIMA(2, 2, 0) RMSE=1346.598
ARIMA(2, 2, 1) RMSE=1039.660
ARIMA(3, 0, 0) RMSE=959.165
ARIMA(3, 1, 0) RMSE=1028.226
ARIMA(3, 1, 1) RMSE=979.940
ARIMA(3, 2, 0) RMSE=1244.846
ARIMA(4, 0, 0) RMSE=968.411
ARIMA(4, 1, 0) RMSE=1046.524
ARIMA(4, 1, 1) RMSE=1006.261
ARIMA(4, 2, 0) RMSE=1227.556
ARIMA(4, 2, 1) RMSE=1055.488
ARIMA(5, 0, 0) RMSE=987.690
ARIMA(5, 1, 0) RMSE=1042.642
ARIMA(5, 1, 1) RMSE=1023.402
ARIMA(5, 1, 2) RMSE=1003.159
ARIMA(5, 2, 1) RMSE=1053.743
ARIMA(6, 0, 0) RMSE=996.466
ARIMA(6, 1, 0) RMSE=1018.211
ARIMA(6, 1, 1) RMSE=1023.760
Best ARIMA(0, 0, 1) RMSE=939.[...]<jupyter_text>It seems that the best ARIMA parameters for this problem is (p,d,q)=(0,0,1).In order to review the residual forecast errors of our model, we're going to look at summary statistics of the residuals, a histogram as well as a density plot. We can use the information to modify our model and reduce the error.<jupyter_code>predictions = list()
history = [x for x in train]
for i in range(len(test)):
    diff = difference(history, months_in_year)
    # predict
    model = ARIMA(diff, order=(0,0,1))
    model_fit = model.fit(trend='nc', disp=0)
    yhat = model_fit.forecast()[0]
    yhat = inverse_difference(history, yhat, months_in_year)
    predictions.append(yhat)
    # observation
    obs = test[i]
    history.append(obs)
    
# residual errors
residuals = [test[i]-predictions[i] for i in range(len(test))]
residuals = pd.DataFrame(residuals)
print(residuals.describe())
# plot
residuals.hist(ax=plt.gca())
plt.show()
residuals.plot(kind='kde', ax=plt.gca())
plt.show()
bias = 165.904730
predictions = list()
history = [x for x in train]
for i in range(len(test)):
    diff = difference(history, months_in_year)
    # predict
    model = ARIMA(diff, order=(0,0,1))
    model_fit = model.fit(trend='nc', disp=0)
    yhat = model_fit.forecast()[0]
    yhat = bias + inverse_difference(history, yhat, months_in_year)
    predictions.append(yhat)
    # observation
    obs = test[i]
    history.append(obs)
    
# residual errors
residuals = [test[i]-predictions[i] for i in range(len(test))]
residuals = pd.DataFrame(residuals)
rmse = sqrt(mean_squared_error(test, predictions))
print('RMSE: %.3f' % (rmse))
print(residuals.describe())
# residual errors
residuals = [test[i]-predictions[i] for i in range(len(test))]
residuals = pd.DataFrame(residuals)
# plot
plot_acf(residuals, ax=plt.gca())
plt.show()
plot_pacf(residuals, ax=plt.gca())
plt.show()<jupyter_output><empty_output><jupyter_text>Model validationWe will now train our ARIMA model one last time and save it in order to test its performance on the validation set (the final 12 months) we seperated from the original dataset.Afterwards we'll plot our model's prediction with the actual sales to see visually how well our model performs.<jupyter_code>from scipy.stats import boxcox
from statsmodels.tsa.arima_model import ARIMAResults

# monkey patch for the ARIMA class due to bug
def __getnewargs__(self):
    return ((self.endog), (self.k_lags, self.k_diff, self.k_ma))

ARIMA.__getnewargs__=__getnewargs__

# load data
series = pd.read_csv('dataset.csv', index_col=0)
X = series.values.astype('float32')
diff = difference(X, months_in_year)

# fit model
model = ARIMA(diff, order=(0,0,1))
model_fit = model.fit(trend='nc', disp=0)

# bias (from residual mean)
bias = 165.904730

# save model
model_fit.save('model.pkl')
np.save('model_bias.npy', [bias])

# prepare datasets
history = [x for x in X]
validation = pd.read_csv('validation.csv', index_col=0)
y = validation.values.astype('float32')

# we make the first prediction to begin the forecasting
# and then we retrain our model with the new timestep. 
# We loop this process until we reach the end of our
# validation set.
predictions = list()
yhat = float(model_fit.forecast()[0])
yhat = bias + inverse_difference(history, yhat, months_in_year)
predictions.append(yhat)
history.append(y[0])
print('Sales: Predicted = %.3f million, Expected = %.3f million' % (yhat, y[0]))
# rolling forecasts
for i in range(1, len(y)):
    diff = difference(history, months_in_year)
    # predict
    model = ARIMA(diff, order=(0,0,1))
    model_fit = model.fit(trend='nc', disp=0)
    yhat = model_fit.forecast()[0]
    yhat = bias + inverse_difference(history, yhat, months_in_year)
    predictions.append(yhat)
    # observation
    obs = y[i]
    history.append(obs)
    print('Sales: Predicted = %.3f million, Expected = %.3f million' % (yhat, y[i]))
    
rmse = sqrt(mean_squared_error(y, predictions))
print('RMSE: %.3f' % rmse)
plt.plot(y)
plt.plot(predictions, color='red')
plt.show()<jupyter_output>Sales: Predicted = 6794.773 million, Expected = 6981.000 million
Sales: Predicted = 10101.763 million, Expected = 9851.000 million
Sales: Predicted = 13219.067 million, Expected = 12670.000 million
Sales: Predicted = 3996.536 million, Expected = 4348.000 million
Sales: Predicted = 3465.934 million, Expected = 3564.000 million
Sales: Predicted = 4522.684 million, Expected = 4577.000 million
Sales: Predicted = 4901.336 million, Expected = 4788.000 million
Sales: Predicted = 5190.094 million, Expected = 4618.000 million
Sales: Predicted = 4930.190 million, Expected = 5312.000 million
Sales: Predicted = 4944.785 million, Expected = 4298.000 million
Sales: Predicted = 1699.409 million, Expected = 1413.000 million
Sales: Predicted = 6085.324 million, Expected = 5877.000 million
RMSE: 361.110<jupyter_text>Looks like ARIMA works very well at predicting the champagne sales.  MLP ModelWe will now use a multilayer perceptron in order to compare it to ARIMA. We've already created a differenced series in order to make the time series stationary. In the case of the MLP model, we'll be also turning the problem into a supervised learning problem and rescale the data to values between -1 and 1.<jupyter_code>from sklearn.preprocessing import MinMaxScaler
from keras.models import Sequential
from keras.layers import Dense
import matplotlib
# be able to save images on server
matplotlib.use('Agg')
from matplotlib import pyplot as plt

def scale(train, test):
    # fit scaler
    scaler = MinMaxScaler(feature_range=(-1, 1))
    scaler = scaler.fit(train)
    # transform train
    train = train.reshape(train.shape[0], train.shape[1])
    train_scaled = scaler.transform(train)
    # transform test
    test = train.reshape(test.shape[0], test.shape[1])
    test_scaled = scaler.transform(test)
    return scaler, train_scaled, test_scaled

# invert scaling for a forecasted value (after we've made prediction)
def invert_scale(scaler, X, yhat):
    new_row = [x for x in X] + [yhat]
    array = numpy.array(new_row)
    array = array.reshape(1, len(array))
    inverted = scaled.inverse_transform(array)
    return inverted[0, -1]

# MLP model with Keras
def fit_MLP_model(train, batch_size, nb_epoch, neurons):
    # Create model
    X, y = train[:, 0:-1], train[:, -1]
    model = Sequential()
    model.add(Dense(neurons, activation='relu', input_dim=X.shape[1]))
    model.add(Dense(1))
    # Compile model
    model.compile(loss='mean_squared_error', optimizer='adam')
    model.fit(X, y, epochs=nb_epochs, batch_size=batch_size, verbose=0, shuffle=False)
    return model<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11778162482997404, 'fraction_numerical': 0.08507481142574502, 'mean_word_length': 4.020173805090006, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 37, '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': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '572072', 'license': 'MIT', 'n_tokens_mistral': 6716, 'n_tokens_neox': 5742, 'n_words': 1625, 'path': 'Monthly Champagne Sales/Monthly-Champagne-Sales.ipynb', 'repo_name': 'JayThibs/time-series-forecasting', 'size': '16174'}