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
