---

| Top | SKTime Overview | Traffic Dataset | Pre-processing | FB Prophet Model | Forecasting Function | AutoARIMA Model |

---

Forecasting with SKTime (with Dave Ebbelaar)
| YouTube Tutorial | Tutorials GitHub | Dave Ebbelaar |

Evan Marie online: | EvanMarie.com| EvanMarie@Proton.me | Linked In | GitHub | Hugging Face | Mastadon | Jovian.ai | TikTok | CodeWars | Discord ⇨ ✨ EvanMarie ✨#6114 |

from helpers_strawberries import *
import_all()
plt.style.use('strawberries.mplstyle')
%matplotlib inline

%%html
<style>
a:link {color: #BF4141 !important; font-weight: 600 !important;}
a:visited {color: #BF4141 !important; font-weight: 600 !important;}
</style>

---

| Top | SKTime Overview | Traffic Dataset | Pre-processing | FB Prophet Model | Forecasting Function | AutoARIMA Model |

---

SKTime Overview
- SKTime example - Using the airline dataset from SKTime - Running a forecasting model and predicting

from sktime.datasets import load_airline
from sktime.forecasting.base import ForecastingHorizon as FH
from sktime.forecasting.model_selection import temporal_train_test_split as split
from sktime.forecasting.theta import ThetaForecaster
from sktime.performance_metrics.forecasting import mean_absolute_percentage_error as mape
from sktime.utils.plotting import plot_series
import seaborn as sns

# investigating imported data

data = load_airline()
data.plot(title = 'SKTime Airline Time Series');

temporal_train_test_split()- SKTime splitter makes sure training data always comes before testing

training, testing = split(data, train_size = .75)
pretty(len(training), 'len(training)')
pretty(len(testing), 'len(testing)')

len(training)
108

len(testing)
36

ForecastingHorizon() - this will match the indices of the test data, that which we want to predict

horizon = FH(testing.index, is_relative = False)

pretty(horizon, 'horizon = FH(testing.index, is_relative = False)')

horizon = FH(testing.index, is_relative = False)

ForecastingHorizon(['1958-01', '1958-02', '1958-03', '1958-04', '1958-05', '1958-06', '1958-07', '1958-08', '1958-09', '1958-10', '1958-11', '1958-12', '1959-01', '1959-02', '1959-03', '1959-04', '1959-05', '1959-06', '1959-07', '1959-08', '1959-09', '1959-10', '1959-11', '1959-12', '1960-01', '1960-02', '1960-03', '1960-04', '1960-05', '1960-06', '1960-07', '1960-08', '1960-09', '1960-10', '1960-11', '1960-12'], dtype='period[M]', name='Period', is_relative=False)

ThetaForecaster() - Season period is set to 12, which is 12 months, 1 year - monthly, seasonly, periodically - 1 indicates yearly, 12 indicates monthly, and 4 indicates quarterly

caster = ThetaForecaster(sp = 12)
caster.fit(training)

ThetaForecaster(sp=12)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

# getting predictions

predictions = caster.predict(horizon)
mape_score = mape(testing, predictions)
pretty(f'{mape_score * 100: .3f}%', 
	   'Mean Absolute Precentage Error', fontsize = '13px')

Mean Absolute Precentage Error
8.175%

# predictions versus actual numbers

plot_series(data, predictions, 
			colors = ['#A3D9C9','#BF656D'], 
		   title = 'Training, Testing, and Predictions');

---

| Top | SKTime Overview | Traffic Dataset | Pre-processing | FB Prophet Model | Forecasting Function | AutoARIMA Model |

---

Traffic Dataset
- importing the traffic dataset - Dataset on Kaggle

data = pd.read_csv('https://raw.githubusercontent.com/daveebbelaar/data-science-tutorials/main/data/raw/traffic.csv', parse_dates = [0], index_col = [0])

# dataset overview

df_overview(data, 'Traffic Dataset', 4, )

Traffic Dataset Columns

	Junction	Vehicles	ID
datatype	int64	int64	int64
missing values	0	0	0
count	48,120.00	48,120.00	48,120.00
mean	2.18	22.79	20,163,302,501.22
std	0.97	20.75	5,944,853.80
min	1.00	1.00	20,151,101,001.00
25%	1.00	9.00	20,160,416,019.00
50%	2.00	15.00	20,160,930,037.00
75%	3.00	29.00	20,170,225,162.25
max	4.00	180.00	20,170,630,234.00

   

Traffic Dataset Key Points

total rows	48,120
total columns	3
column names	Junction, Vehicles, ID
index start	2015-11-01 00:00:00
index end	2017-06-30 23:00:00
total missing values	0

   

Traffic Dataset Head and Tail

head(4)

	Junction	Vehicles	ID
DateTime
2015-11-01	1	15	20,151,101,001
2015-11-01	1	13	20,151,101,011
2015-11-01	1	10	20,151,101,021
2015-11-01	1	7	20,151,101,031

   

tail(4)

	Junction	Vehicles	ID
DateTime
2017-06-30	4	30	20,170,630,204
2017-06-30	4	16	20,170,630,214
2017-06-30	4	22	20,170,630,224
2017-06-30	4	12	20,170,630,234

   

# pivot() - so that the columns are the 4 junctions, and the values are the numbers of vehicles at each in any given hour

data = data.pivot(columns = 'Junction', values = 'Vehicles')
head_tail_horz(data, 5, 'Traffic Data Pivoted')

Traffic Data Pivoted

head(5)

Junction	1	2	3	4
DateTime
2015-11-01	15.00	6.00	9.00	nan
2015-11-01	13.00	6.00	7.00	nan
2015-11-01	10.00	5.00	5.00	nan
2015-11-01	7.00	6.00	1.00	nan
2015-11-01	9.00	7.00	2.00	nan

   

tail(5)

Junction	1	2	3	4
DateTime
2017-06-30	105.00	34.00	33.00	11.00
2017-06-30	96.00	35.00	31.00	30.00
2017-06-30	90.00	31.00	28.00	16.00
2017-06-30	84.00	29.00	26.00	22.00
2017-06-30	78.00	27.00	39.00	12.00

   

# visualization - the hourly traffic at the 4 different junctions

data.plot(title = 'Hourly Traffic 2016-2017');

---

| Top | SKTime Overview | Traffic Dataset | Pre-processing | FB Prophet Model | Forecasting Function | AutoARIMA Model |

---

Pre-processing the Data
- The data requires resampling, because Pandas does not recognize the frequency of the timestamps

resample() for frequency - Pandas does not recognize the frequency or timestep between each of the index entries - rule = 'H' - resamples to hourly - sum() - we want to know the total number of vehicles per hour

data.index

DatetimeIndex(['2015-11-01 00:00:00', '2015-11-01 01:00:00',
               '2015-11-01 02:00:00', '2015-11-01 03:00:00',
               '2015-11-01 04:00:00', '2015-11-01 05:00:00',
               '2015-11-01 06:00:00', '2015-11-01 07:00:00',
               '2015-11-01 08:00:00', '2015-11-01 09:00:00',
               ...
               '2017-06-30 14:00:00', '2017-06-30 15:00:00',
               '2017-06-30 16:00:00', '2017-06-30 17:00:00',
               '2017-06-30 18:00:00', '2017-06-30 19:00:00',
               '2017-06-30 20:00:00', '2017-06-30 21:00:00',
               '2017-06-30 22:00:00', '2017-06-30 23:00:00'],
              dtype='datetime64[ns]', name='DateTime', length=14592, freq=None)

data.resample(rule = 'H').sum().index

DatetimeIndex(['2015-11-01 00:00:00', '2015-11-01 01:00:00',
               '2015-11-01 02:00:00', '2015-11-01 03:00:00',
               '2015-11-01 04:00:00', '2015-11-01 05:00:00',
               '2015-11-01 06:00:00', '2015-11-01 07:00:00',
               '2015-11-01 08:00:00', '2015-11-01 09:00:00',
               ...
               '2017-06-30 14:00:00', '2017-06-30 15:00:00',
               '2017-06-30 16:00:00', '2017-06-30 17:00:00',
               '2017-06-30 18:00:00', '2017-06-30 19:00:00',
               '2017-06-30 20:00:00', '2017-06-30 21:00:00',
               '2017-06-30 22:00:00', '2017-06-30 23:00:00'],
              dtype='datetime64[ns]', name='DateTime', length=14592, freq='H')

daily - if the resampling frequency is set to daily, the data is aggregated for each column by day - This results in far fewer rows

traffic_daily = data.resample(rule = 'D').sum()
head_tail_horz(traffic_daily, 5, 'Traffic Data with Daily Frequency')

Traffic Data with Daily Frequency

head(5)

Junction	1	2	3	4
DateTime
2015-11-01	327.00	133.00	136.00	0.00
2015-11-02	546.00	197.00	166.00	0.00
2015-11-03	544.00	217.00	150.00	0.00
2015-11-04	498.00	199.00	121.00	0.00
2015-11-05	464.00	200.00	106.00	0.00

   

tail(5)

Junction	1	2	3	4
DateTime
2017-06-26	1,774.00	587.00	414.00	181.00
2017-06-27	2,187.00	730.00	567.00	225.00
2017-06-28	2,080.00	720.00	503.00	276.00
2017-06-29	2,086.00	696.00	497.00	296.00
2017-06-30	1,883.00	655.00	553.00	272.00

   

---

| Top | SKTime Overview | Traffic Dataset | Pre-processing | FB Prophet Model | Forecasting Function | AutoARIMA Model |

---

Forecasting with FB Prophet
"Prophet is a procedure for forecasting time series data based on an additive model where non-linear trends are fit with yearly, weekly, and daily seasonality, plus holiday effects. It works best with time series that have strong seasonal effects and several seasons of historical data. Prophet is robust to missing data and shifts in the trend, and typically handles outliers well." (Source)

from sklearn.metrics import mean_absolute_error, r2_score
from sktime.forecasting.base import ForecastingHorizon
from sktime.forecasting.fbprophet import Prophet

forecaster = Prophet()

horizon = 30 - predicting 30 days into the future

horizon = 30

Working with just one of the four junctions at first

single = traffic_daily[1]
single_training = single[:-horizon]
single_testing = single.tail(horizon)

forecaster.fit(single_training)

15:32:46 - cmdstanpy - INFO - Chain [1] start processing
15:32:46 - cmdstanpy - INFO - Chain [1] done processing

Prophet()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

forecast_horizon = ForecastingHorizon(single_testing.index, is_relative = False)

pretty(forecast_horizon, 'forecast_horizon')

forecast_horizon

ForecastingHorizon(['2017-06-01', '2017-06-02', '2017-06-03', '2017-06-04', '2017-06-05', '2017-06-06', '2017-06-07', '2017-06-08', '2017-06-09', '2017-06-10', '2017-06-11', '2017-06-12', '2017-06-13', '2017-06-14', '2017-06-15', '2017-06-16', '2017-06-17', '2017-06-18', '2017-06-19', '2017-06-20', '2017-06-21', '2017-06-22', '2017-06-23', '2017-06-24', '2017-06-25', '2017-06-26', '2017-06-27', '2017-06-28', '2017-06-29', '2017-06-30'], dtype='datetime64[ns]', name='DateTime', freq='D', is_relative=False)

predictions: number of vehicles as well as confidence intervals

predictions = forecaster.predict(forecast_horizon)
see(predictions.head(5), 'first five predictions',
    fontsize = '13px')

first five predictions

	0
2017-06-01	1,789.96
2017-06-02	1,718.58
2017-06-03	1,430.59
2017-06-04	1,397.08
2017-06-05	1,793.11

confidence interval predictions (90% coverage)
In frequentist statistics, a confidence interval (CI) is a range of estimates for an unknown parameter. A confidence interval is computed at a designated confidence level; the 95% confidence level is most common, but other levels, such as 90% or 99%, are sometimes used. The confidence level represents the long-run proportion of CIs (at the given confidence level) that theoretically contain the true value of the parameter. For example, out of all intervals computed at the 95% level, 95% of them should contain the parameter's true value. Factors affecting the width of the CI include the sample size, the variability in the sample, and the confidence level. All else being the same, a larger sample produces a narrower confidence interval, greater variability in the sample produces a wider confidence interval, and a higher confidence level produces a wider confidence interval.
(Source)

confidence = forecaster.predict_interval(forecast_horizon, coverage = 0.9)
see(confidence.head(5), 'first five confidence interval predictions', 
	fontsize = '13px')

first five confidence interval predictions

	Coverage
	0.90
	lower	upper
2017-06-01	1,581.58	1,965.31
2017-06-02	1,527.12	1,903.89
2017-06-03	1,258.70	1,615.81
2017-06-04	1,204.58	1,579.13
2017-06-05	1,594.32	1,990.08

actuals = single.tail(horizon)
mean_abs_error = mean_absolute_error(actuals, predictions)
pretty(mean_abs_error, 'Mean Absolute Error for Predictions')

Mean Absolute Error for Predictions
184.347541

pretty(f'{mean_abs_error / actuals.mean() * 100: .3f}%',
	   'Relative percentage error')

Relative percentage error
10.464%

plt.plot(single.tail(horizon * 3),
		 label = 'actuals',
		 color = 'cyan');

plt.gca().fill_between(confidence.index, 
					  (confidence.iloc[:, 0]),
					  (confidence.iloc[:, 1]),
					  color = 'white', alpha = 0.3);

plt.plot(predictions, label = 'predictions', color = 'yellow')
plt.ylim(bottom = 1000)
plt.legend(loc = 2);
plt.title('Actuals vs Predictions (with confidence)');

Making Future Predictions
- use all data to train, since this will go into the future - using the same 30 day forecast horizon -

forecaster.fit(single)

15:36:27 - cmdstanpy - INFO - Chain [1] start processing
15:36:27 - cmdstanpy - INFO - Chain [1] done processing

Prophet()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

data_endpoint = single.index.max()
pretty(data_endpoint, 'End of training data:')

End of training data:
2017-06-30 00:00:00

forecast_horizon = ForecastingHorizon(pd.date_range(str(data_endpoint), 
												 periods = horizon, 
												 freq = 'D'),
								   is_relative = False)

Forecasting Horizon - this will go from the end of the data input and one month into the future - this creates the date range

pretty(forecast_horizon, 'Forecasting Horizon')

Forecasting Horizon

ForecastingHorizon(['2017-06-30', '2017-07-01', '2017-07-02', '2017-07-03', '2017-07-04', '2017-07-05', '2017-07-06', '2017-07-07', '2017-07-08', '2017-07-09', '2017-07-10', '2017-07-11', '2017-07-12', '2017-07-13', '2017-07-14', '2017-07-15', '2017-07-16', '2017-07-17', '2017-07-18', '2017-07-19', '2017-07-20', '2017-07-21', '2017-07-22', '2017-07-23', '2017-07-24', '2017-07-25', '2017-07-26', '2017-07-27', '2017-07-28', '2017-07-29'], dtype='datetime64[ns]', freq='D', is_relative=False)

predictions = forecaster.predict(forecast_horizon)
confidence = forecaster.predict_interval(forecast_horizon, 
										coverage = 0.9).astype('float')

plt.plot(single.tail(horizon * 3),
		 label = 'actuals',
		 color = 'cyan');

plt.gca().fill_between(confidence.index, 
					  (confidence.iloc[:, 0]),
					  (confidence.iloc[:, 1]),
					  color = 'white', alpha = 0.3);
plt.plot(predictions, label = 'predictions', color = 'yellow');
plt.ylim(bottom = 800);
plt.xticks(fontsize = 12);
plt.title('Forecasting into Future')

Text(0.5, 1.0, 'Forecasting into Future')

---

| Top | SKTime Overview | Traffic Dataset | Pre-processing | FB Prophet Model | Forecasting Function | AutoARIMA Model |

---

Forecasting Function
- loop over a dataframe of time series data - train an SKTime forecasting model - make predictions - visualize the results

def forecast_data(df, horizon, forecaster, 
				  future = True, 
				  confidence_coverage = 0.9,
				  frequency = 'D'):
	
	forecast_data = df.resample(rule = frequency).sum()
	
	#interpolate any missing data
	forecast_data = forecast_data.interpolate(method = 'time')
	
	for column in df.columns:
		
		if future:
			current = forecast_data[column].dropna()
			forecaster.fit(current)
			data_endpoint = current.index.max()
			forecasting_horizon = ForecastingHorizon(
									pd.date_range(str(data_endpoint),
												 periods = horizon,
												 freq = frequency),
									is_relative = False)
			
			predictions = forecaster.predict(forecasting_horizon)
			confidence = forecaster.predict_interval(forecasting_horizon, 
													 coverage = confidence_coverage)\
													 .astype('float')
			mean_abs_error = np.NaN
			
		
		else:
			current = forecast_data[column]
			training = current[:-horizon]
			testing = current.tail(horizon)
			
			forecaster.fit(training)
			forecasting_horizon = ForecastingHorizon(testing.index, 
													 is_relative = False)
			predictions = forecaster.predict(forecasting_horizon)
			confidence = forecaster.predict_interval(forecasting_horizon,
													 coverage = confidence_coverage)\
													 .astype('float')
			actuals = current.tail(horizon)
			mean_abs_error = mean_absolute_error(actuals, predictions)
			
		# Visualize results
		plt.plot(current.tail(horizon * 3), label = 'actual', color = 'cyan')
		plt.gca().fill_between(confidence.index, 
							   (confidence.iloc[:, 0]), (confidence.iloc[:, 1]), 
							   color="white", alpha=0.3)
		plt.plot(predictions, label="predictions", color = 'yellow')
		plt.title(f"{horizon} day forecast for {column} (mae: {round(mean_abs_error, 2)}, confidence: {confidence_coverage*100}%)", size = 22)
		plt.ylim(bottom=100)
		plt.legend(loc = 4)
		plt.grid(False)
		plt.show()

Forecaster - yearly_seasonality and weekly_seasonality are True, because there is an assumption with traffic data that there will be such patterns

forecaster = Prophet(yearly_seasonality = True, weekly_seasonality = True)

forecast_data(traffic_daily, horizon = 30, forecaster = forecaster, future = False)

15:36:42 - cmdstanpy - INFO - Chain [1] start processing
15:36:43 - cmdstanpy - INFO - Chain [1] done processing

15:36:43 - cmdstanpy - INFO - Chain [1] start processing
15:36:43 - cmdstanpy - INFO - Chain [1] done processing

15:36:43 - cmdstanpy - INFO - Chain [1] start processing
15:36:43 - cmdstanpy - INFO - Chain [1] done processing

15:36:43 - cmdstanpy - INFO - Chain [1] start processing
15:36:43 - cmdstanpy - INFO - Chain [1] done processing

---

| Top | SKTime Overview | Traffic Dataset | Pre-processing | FB Prophet Model | Forecasting Function | AutoARIMA Model |

---

Forecasting with AutoARIMA

"An autoregressive integrated moving average, or ARIMA, is a statistical analysis model that uses time-series data to better understand the data set or predict future trends. A statistical model is autoregressive if it predicts future values based on past values. For example, an ARIMA model might seek to predict a stock’s future prices based on its past performance or forecast a company’s earnings based on past periods." (Source)

from sktime.forecasting.arima import AutoARIMA

forecaster = AutoARIMA(sp = 1, suppress_warnings = True)
forecast_data(traffic_daily, horizon = 30, forecaster = forecaster, future = False)

---

| Top | SKTime Overview | Traffic Dataset | Pre-processing | FB Prophet Model | Forecasting Function | AutoARIMA Model |

---