A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

1 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

2 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mlp

#Loading data
data = pd.read_csv('Datasets/AirPassengers.csv',
index_col='Month', parse_dates=['Month'])
data = data.rename(columns
={'#Passengers':'no_passengers'})
data

read_csv().

parse_dates

to_datetime().

index_col='Month'

3 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

plt.plot(data)
plt.show()

4 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

data_for_dist_fitting = data[-70:]

5 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

data_train =
data[~data.isin(data_for_dist_fitting).all(1)]

test_data = data_for_dist_fitting[-20:]

data_for_dist_fitting=data_for_dist_fitting[~data_for_dist
_fitting.isin(test_data).all(1)]

train = plt.plot(data_train,color='blue', label = 'Train

data')

data_f_mc = plt.plot(data_for_dist_fitting, color ='red',

label ='Data for distribution fitting')

test = plt.plot(test_data, color ='black', label = 'Test

data')

plt.legend(loc='best')
plt.title('Data division')
plt.show(block=False)

6 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

7 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

from statsmodels.tsa.stattools import adfuller

def test_stationarity(timeseries):
#Determining rolling statistics
rolmean = timeseries.rolling(window=12).mean()
rolstd = timeseries.rolling(window=12).std()

#plot rolling statistics:

orig = plt.plot(timeseries,color='blue', label =
'Original')
mean = plt.plot(rolmean, color ='red', label ='Rolling
Mean')
std = plt.plot(rolstd, color ='black', label =
'Rolling Std')
plt.legend(loc='best')
plt.title('Rolling Mean and Standard Deviation')
plt.show(block=False)

8 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

#Perform Dickey Fuller test:

print('Results of Dickey Fuller Test:')
dftest = adfuller(timeseries, 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['Critical Value (%s)'%key] = value
print (dfoutput)

test_stationarity(data_train)

9 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

statsmodels seasonal_decompose

from statsmodels.tsa.seasonal import seasonal_decompose

decomposition = seasonal_decompose(data_train)

trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid

plt.subplot(411)
plt.plot(ts_log, label='Original')
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()

10 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

SARIMA(p,d,q).

(P,D,Q).m

11 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

import itertools
p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in
list(itertools.product(p, d, q))]
print('Examples of parameter for SARIMA...')
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))

12 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

import warnings
warnings.filterwarnings('ignore')
import statsmodels.api as sm
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
mod =
sm.tsa.statespace.SARIMAX(data_train,order=param,seasonal_
order=param_seasonal,enforce_stationarity=False,enforce_in
vertibility=False)
results = mod.fit()
print('SARIMA{}x{}12 -
AIC:{}'.format(param,param_seasonal,results.aic))
except Exception as E:
print(E)
continue

13 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

14 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

SARIMAX(1, 1, 1)x(1, 1, 1, 12).

from statsmodels.tsa.statespace.sarimax import SARIMAX

mod= SARIMAX(data_train,order=(1,1,1),seasonal_order=(1,
1, 1, 12),enforce_invertibility=False,
enforce_stationarity=False)
results = mod.fit(disp=0)
print(results.summary())

15 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

pred_sarima = results.forecast(50)
predicted =plt.plot(pred_sarima,label='Prediction by
SARIMA', color='red')
Actual = plt.plot(data_for_dist_fitting,label='Actual
data')
plt.legend(loc='best')
plt.title('SARIMA MODEL')
plt.show(block=False)

16 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

# plot residual errors of the training data

residual_error = pd.DataFrame(results.resid)
residual_error.plot()
plt.show()

residual_error.plot(kind='kde')
plt.show()
print(residual_error.describe())

17 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

#to suppress warnings

warnings.filterwarnings('ignore')

from sklearn.metrics import mean_squared_error

#creating new dataframe for rolling forescast

18 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

history = np.log(data_train.astype(float))
predictions = list()

for i in range(len(data_for_dist_fitting)):
model = SARIMAX(history,order=
(1,1,1),seasonal_order=(1, 1, 1,
12),enforce_invertibility=False,
enforce_stationarity=False)
model_fit = model.fit(disp = 0)
# generate forcecast for next period
output = model_fit.forecast()
#Save the prediction value in yhat
yhat = np.e ** output[0]
#Append yhat to the list of prediction
predictions.append(yhat)
# grabs the observation at the ith index
obs = data_for_dist_fitting[i : i + 1]
# appends the observation to the estimation data set
history = history.append(np.log(obs.astype(float)))

# prints the MSE of the model for the rolling forecast

period
error = mean_squared_error(data_for_dist_fitting,
predictions)
print('Test MSE: %.3f' % error)

# converts the predictions list to a pandas dataframe with

the same index as the actual values
# for plotting purposes
predictions = pd.DataFrame(predictions)
predictions.index = data_for_dist_fitting.index

# sets the plot size to 12x8

mlp.rcParams['figure.figsize'] = (12,8)

# plots the predicted and actual stock prices

plt.plot(data_for_dist_fitting,label='Actual values')
plt.plot(predictions, color = 'red', label='predicted
rolling forecast')
plt.legend(loc='best')
plt.xlabel('week')
plt.ylabel('#passengers')
plt.title('Predicted vs. Actual #of passengers')
plt.show()

19 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

# to suppress warnings
warnings.filterwarnings('ignore')

# sets the plot size to 12x8

mlp.rcParams['figure.figsize'] = (12,8)

# plots the rolling forecast error

rf_errors = data_for_dist_fitting.no_passengers -
predictions[0]
rf_errors.plot(kind = 'kde')

# produces a summary of rolling forecast error

rf_errors.astype(float).describe()

20 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

# to suppress warnings

21 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

import warnings
warnings.filterwarnings('ignore')

# imports the fitter function and produces estimated fits

for our rsarima_errors
from fitter import
Fitter,get_common_distributions,get_distributions

f = Fitter(rf_errors, distributions=
['binomial','norm','laplace','uniform'])
f.fit()
f.summary()

22 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

data_for_dist_fitting

np.random.laplace(loc,scale,size)

def lapace_mc_randv_distribution(mean, rf_errors, n_sim):

#gets the estimated beta or mean absolute distance from

the mean
var = (sum(abs(rf_errors - np.mean(rf_errors)))
/ len(rf_errors))

# uses the numpy function to generate an array of

simulated values
est_range = np.random.laplace(mean,var,n_sim)
# converts the array to a list
est_range = list(est_range)
# returns the simulated values
return(est_range)

23 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

def rolling_forecast_MC(train, test, std_dev, n_sims):

# create a new dataframe that will be added to as the

forecast rolls
history = np.log(data_train.astype(float))

# create an empty list that will hold predictions

predictions = list()

# loops through the indexes of the set being forecasted

for i in range(len(test_data)):
model = SARIMAX(history,order=
(1,1,1),seasonal_order=(1, 1, 1,
12),enforce_invertibility=False,
enforce_stationarity=False)
model_fit = model.fit(disp = 0)

# generate forcecast for next period

output = model_fit.forecast().values

#Save the prediction value in yhat

yhat = np.e ** output[0]

# performs monte carlo simulation using the

predicted price as the mean, user-specified
# standard deviation, and number of simulations
randv_range =
lapace_mc_randv_distribution(yhat,std_dev,n_sims)

#Append yhat to the list of prediction

predictions.append([float(i) for i in
randv_range])

# grabs the observation at the ith index

obs = test_data[i : i + 1]

# appends the observation to the estimation data

set
history =
history.append(np.log(obs.astype(float)))

24 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

# converts the predictions list to a pandas

dataframe with the same index as the actual
# values for plotting purposes
predictions = pd.DataFrame(predictions)
predictions.index = test_data.index

# returns predictions
return(predictions)

data_train = data_train.append(data_for_dist_fitting)

test_preds = rolling_forecast_MC(data_train,
test_data,
rf_errors,
1000)

MC = plt.plot(test_preds)
Actual=plt.plot(test_data,color='black',label='Actual
Demand')
plt.legend(loc='best')
plt.show()

25 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

print('Expected demand:',np.mean(test_preds.values))
print('Quantile(5%):',np.percentile(test_preds,5))
print('Quantile(95%):',np.percentile(test_preds,95))

26 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

predictions.values.tolist().

def rolling_forecast_MC_for_minmax_range(train, test,

std_dev, n_sims):
# create a new dataframe that will be added to as the
forecast rolls
history = np.log(data_train.astype(float))
# create an empty list that will hold predictions
predictions = list()

# loops through the indexes of the set being forecasted

for i in range(len(test_data)):
model = SARIMAX(history,order=
(1,1,1),seasonal_order=(1, 1, 1,
12),enforce_invertibility=False,
enforce_stationarity=False)
model_fit = model.fit(disp = 0)

# generate forcecast for next period

output = model_fit.forecast().values

#Save the prediction value in yhat

yhat = np.e ** output[0]

# performs monte carlo simulation using the

predicted price as the mean, user-specified
# standard deviation, and number of simulations
randv_range =
lapace_mc_randv_distribution(yhat,std_dev,n_sims)

#Append yhat to the list of prediction

predictions.append([float(i) for i in
randv_range])

# grabs the observation at the ith index

obs = test_data[i : i + 1]

# appends the observation to the estimation data set

history =
history.append(np.log(obs.astype(float)))

# converts the predictions list to a pandas dataframe

27 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

with the same index as the actual

# values for plotting purposes

predictions = pd.DataFrame(predictions)
# converts all the estimated yhats in each column to
one list per row
predictions['predicted_range'] =
predictions.values.tolist()

# grabs only the column with all values in a list

predictions =
pd.DataFrame(predictions['predicted_range'])
predictions.index = test_data.index

# returns predictions
return(predictions)

# produces a rolling forecast with prediction intervals

using 1000 MC sims
test_preds_minmax =
rolling_forecast_MC_for_minmax_range(data_train,
test_data,
rf_errors,
1000)

test_preds_minmax.head()

28 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

# creates an empty list

prediction_interval = []
# loops through the rows in the testing data set
for i in range(len(test_data)):
# appends true if the actual price is in the interval
of predicted prices and false
# otherwise

prediction_interval.append(np.where(min(test_preds_minmax.
predicted_range[i]) <=

test_data.no_passengers[i]
<=
max(test_preds_minmax.predicted_range[i]),
True, False))
# prints the percentage of actual prices in the prediction
intervals
print('Percentage of Demand in Predicted Demand Range: %f'
%
(100 * sum(prediction_interval) /
len(prediction_interval)))

29 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

# creates empty lists to append to with minimum and maximum

values for each weeks prediction
min_range = []
max_range = []

# loops through the rows in test_preds

for i in range(len(test_preds_minmax)):
# appends to the list the min or max value as appropriate
min_range.append(min(test_preds_minmax.predicted_range[i]))
max_range.append(max(test_preds_minmax.predicted_range[i]))

# converts the lists to data frames and makes their indexes

match up with the dates they're
# predicting
min_range = pd.DataFrame(min_range)
min_range.index = test_data.index
max_range = pd.DataFrame(max_range)
max_range.index = test_data.index

# plots the actual stock price with prediction intervals

plt.plot(test_data, color ='red',label='Actual Data')
plt.plot(min_range, color = 'm', label='Min range')
plt.plot(max_range, color = 'b', label ='Max range')
plt.legend(loc='best')
plt.xlabel('Month')
plt.ylabel('No of Passengers')
plt.title('Actual Demand with Prediction Intervals')
plt.show()

30 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

31 of 32 3/15/2022, 9:49 AM
A Stochastic Model For Demand Forecating In Python | by Ju... https://medium.com/mlearning-ai/a-stochastic-model-for-de...

32 of 32 3/15/2022, 9:49 AM