#Introduction to Time Series Forecasting with Python
By Prof. Angel Colmenares
The goal of this codes is to show you how to get results on univariate time series forecasting problems using the Python ecosystem. This goal cannot be achieved until you apply the lessons from this book on your own projects and get results. This is a guidebook or a cookbook designed for immediate use
Quarto enables you to weave together content and executable code into a finished document. To learn more about Quarto see https://quarto.org.
Chapter I traditional models on Python
Index
1.1 Basics concepts:
Create a rolling mean feature, statistics feature
Create expanding window features
Create a line plot, dot plot, stacked line plots, histogram plot, density plot
1.2 Plots of:
Box and Whisker Plots
Heat Maps
Lag Scatter Plots
Autocorrelation Plots
1.3 Resampling and Interpolation, Upsampling and Downsampling Data
1.4 Power Transforms, Square Root Transform, Log Transform, Box-Cox Transform
1.5 Moving Average Smoothing , Introduction to White Noise, Random Walk
1.6 Decompose Time Series Data
1.7 Use and Remove Trends, Seasonality , Stationarity in Time Series
1.8 Data Backtest Forecast Models, Forecasting Performance Measures
1.9 Persistence Model for Forecasting
1.10 Autoregression Models for Forecasting
1.11 Visualize Residual Forecast Errors
Chapter II Time Series Forecasting Using Generative AI: Leveraging AI for Precision Forecasting on Python
2.1 Multilayer Perceptron (MLP)
2.2 Temporal Convolutional Networks (TCN)
2.3 Bidirectional Temporal Convolutional Network (BiTCN)
2.4 Recurrent Neural Network (RNN)
2.5 Long Short-Term Memory (LSTM)
2.6 Neural Networks Based on Autoregression (DeepAR)
2.7 Neural Basis Expansion Analysis (NBEATS)
2.8 Transformers for Time Series
Vanilla Transformer
Inverted Transformers
2.9 DLinear
2.10 NLinear
2.11 Patch Time Series Transformer (PatchTST)
2.12 Reprogramming Large Language Model (Time-LLM)
Basics concepts:
#!pip install
from pandas import read_csv
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
print(type(series))<class 'pandas.core.series.Series'>print(series.head(10))Date
1959-01-01 35
1959-01-02 32
1959-01-03 30
1959-01-04 31
1959-01-05 44
1959-01-06 29
1959-01-07 45
1959-01-08 43
1959-01-09 38
1959-01-10 27
Name: Births, dtype: int64print(series.size)365print(series['1959-01'])Date
1959-01-01 35
1959-01-02 32
1959-01-03 30
1959-01-04 31
1959-01-05 44
1959-01-06 29
1959-01-07 45
1959-01-08 43
1959-01-09 38
1959-01-10 27
1959-01-11 38
1959-01-12 33
1959-01-13 55
1959-01-14 47
1959-01-15 45
1959-01-16 37
1959-01-17 50
1959-01-18 43
1959-01-19 41
1959-01-20 52
1959-01-21 34
1959-01-22 53
1959-01-23 39
1959-01-24 32
1959-01-25 37
1959-01-26 43
1959-01-27 39
1959-01-28 35
1959-01-29 44
1959-01-30 38
1959-01-31 24
Name: Births, dtype: int64print(series.describe())count 365.000000
mean 41.980822
std 7.348257
min 23.000000
25% 37.000000
50% 42.000000
75% 46.000000
max 73.000000
Name: Births, dtype: float64from pandas import read_csv
from pandas import DataFrame
from pandas import concat
series = read_csv('daily-minimum-temperatures.csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
dataframe = DataFrame()
dataframe['month'] = [series.index[i].month for i in range(len(series))]
dataframe['day'] = [series.index[i].day for i in range(len(series))]
dataframe['temperature'] = [series[i] for i in range(len(series))]<string>:1: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`print(dataframe.head(5)) month day temperature
0 1 1 20.7
1 1 2 17.9
2 1 3 18.8
3 1 4 14.6
4 1 5 15.8temps = DataFrame(series.values)
dataframe = concat([temps.shift(3), temps.shift(2), temps.shift(1), temps], axis=1)
dataframe.columns = ['t-2', 't-1', 't', 't+1']
print(dataframe.head(5)) t-2 t-1 t t+1
0 NaN NaN NaN 20.7
1 NaN NaN 20.7 17.9
2 NaN 20.7 17.9 18.8
3 20.7 17.9 18.8 14.6
4 17.9 18.8 14.6 15.8# create a rolling mean feature
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
series = read_csv('daily-minimum-temperatures.csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
temps = DataFrame(series.values)
shifted = temps.shift(1)
window = shifted.rolling(window=2)
means = window.mean()
dataframe = concat([means, temps], axis=1)
dataframe.columns = ['mean(t-1,t)', 't+1']
print(dataframe.head(5)) mean(t-1,t) t+1
0 NaN 20.7
1 NaN 17.9
2 19.30 18.8
3 18.35 14.6
4 16.70 15.8temps = DataFrame(series.values)
# create rolling statistics feature
width = 3
shifted = temps.shift(width - 1)
window = shifted.rolling(window=width)
dataframe = concat([window.min(), window.mean(), window.max(), temps], axis=1)
dataframe.columns = ['min', 'mean', 'max', 't+1']
print(dataframe.head(5)) min mean max t+1
0 NaN NaN NaN 20.7
1 NaN NaN NaN 17.9
2 NaN NaN NaN 18.8
3 NaN NaN NaN 14.6
4 17.9 19.133333 20.7 15.8# create expanding window features
window = temps.expanding()
dataframe = concat([window.min(), window.mean(), window.max(), temps.shift(-1)], axis=1)
dataframe.columns = ['min', 'mean', 'max', 't+1']
print(dataframe.head(5)) min mean max t+1
0 20.7 20.700000 20.7 17.9
1 17.9 19.300000 20.7 18.8
2 17.9 19.133333 20.7 14.6
3 14.6 18.000000 20.7 15.8
4 14.6 17.560000 20.7 15.8# create a line plot
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
series.plot()
pyplot.show()
# create a dot plot
series.plot(style='k.')
pyplot.show()
# create stacked line plots
from pandas import read_csv
from pandas import DataFrame
from pandas import Grouper
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
groups = series.groupby(Grouper(freq='Y'))
years = DataFrame()
for name, group in groups:
years[name.year] = group.values
years.plot(subplots=True, legend=False)array([<Axes: >, <Axes: >, <Axes: >, <Axes: >, <Axes: >, <Axes: >,
<Axes: >, <Axes: >, <Axes: >, <Axes: >], dtype=object)pyplot.show()
pyplot.close()
# create a histogram plot
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
series.hist()
pyplot.show()
pyplot.close()
# create a density plot
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
series.plot(kind='kde')
pyplot.show()
pyplot.close()
## 6.5 Box and Whisker Plots by Interval
# create a boxplot of yearly data
from pandas import read_csv
from pandas import DataFrame
from pandas import Grouper
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
groups = series.groupby(Grouper(freq='Y'))
years = DataFrame()
for name, group in groups:
years[name.year] = group.values
years.boxplot()
pyplot.show()
pyplot.close()
# create a boxplot of monthly data
from pandas import read_csv
from pandas import DataFrame
from pandas import Grouper
from matplotlib import pyplot
from pandas import concat
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
one_year = series['1990']
groups = one_year.groupby(Grouper(freq='M'))
months = concat([DataFrame(x[1].values) for x in groups], axis=1)
months = DataFrame(months)
months.columns = range(1,13)
months.boxplot()
pyplot.show()
pyplot.close()
## 6.6 Heat Maps
# create a heat map of yearly data
from pandas import read_csv
from pandas import DataFrame
from pandas import Grouper
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
groups = series.groupby(Grouper(freq='Y'))
years = DataFrame()
for name, group in groups:
years[name.year] = group.values
years = years.T
pyplot.matshow(years, interpolation=None, aspect='auto')
pyplot.show()
pyplot.close()
# create a heat map of monthly data
from pandas import read_csv
from pandas import DataFrame
from pandas import Grouper
from matplotlib import pyplot
from pandas import concat
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
one_year = series['1990']
groups = one_year.groupby(Grouper(freq='M'))
months = concat([DataFrame(x[1].values) for x in groups], axis=1)
months = DataFrame(months)
months.columns = range(1,13)
pyplot.matshow(months, interpolation=None, aspect='auto')
pyplot.show()
pyplot.close()
## 6.7 Lag Scatter Plots
# create a scatter plot
from pandas import read_csv
from matplotlib import pyplot
from pandas.plotting import lag_plot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
lag_plot(series)
pyplot.show()
# create multiple scatter plots
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
values = DataFrame(series.values)
lags = 7
columns = [values]
for i in range(1,(lags + 1)):
columns.append(values.shift(i))
dataframe = concat(columns, axis=1)
columns = ['t']
for i in range(1,(lags + 1)):
columns.append('t-' + str(i))
dataframe.columns = columns
pyplot.figure(1)
for i in range(1,(lags + 1)):
ax = pyplot.subplot(240 + i)
ax.set_title('t vs t-' + str(i))
pyplot.scatter(x=dataframe['t'].values, y=dataframe['t-'+str(i)].values)
pyplot.show()
## 6.8 Autocorrelation Plots
# create an autocorrelation plot
from pandas import read_csv
from matplotlib import pyplot
from pandas.plotting import autocorrelation_plot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
autocorrelation_plot(series)
pyplot.show()
pyplot.close()### Chapter 7 Resampling and Interpolation
## 7.3 Upsampling Data
# upsample to daily intervals
from pandas import read_csv
from datetime import datetime
def parser(x):
return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")<string>:1: FutureWarning: The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.upsampled = series.resample('D').mean()
print(upsampled.head(32))Month
1901-01-01 266.0
1901-01-02 NaN
1901-01-03 NaN
1901-01-04 NaN
1901-01-05 NaN
1901-01-06 NaN
1901-01-07 NaN
1901-01-08 NaN
1901-01-09 NaN
1901-01-10 NaN
1901-01-11 NaN
1901-01-12 NaN
1901-01-13 NaN
1901-01-14 NaN
1901-01-15 NaN
1901-01-16 NaN
1901-01-17 NaN
1901-01-18 NaN
1901-01-19 NaN
1901-01-20 NaN
1901-01-21 NaN
1901-01-22 NaN
1901-01-23 NaN
1901-01-24 NaN
1901-01-25 NaN
1901-01-26 NaN
1901-01-27 NaN
1901-01-28 NaN
1901-01-29 NaN
1901-01-30 NaN
1901-01-31 NaN
1901-02-01 145.9
Freq: D, Name: Sales, dtype: float64# upsample to daily intervals with linear interpolation
from pandas import read_csv
from datetime import datetime
from matplotlib import pyplot
def parser(x):
return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")<string>:1: FutureWarning: The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.upsampled = series.resample('D').mean()
interpolated = upsampled.interpolate(method='linear')
print(interpolated.head(32))Month
1901-01-01 266.000000
1901-01-02 262.125806
1901-01-03 258.251613
1901-01-04 254.377419
1901-01-05 250.503226
1901-01-06 246.629032
1901-01-07 242.754839
1901-01-08 238.880645
1901-01-09 235.006452
1901-01-10 231.132258
1901-01-11 227.258065
1901-01-12 223.383871
1901-01-13 219.509677
1901-01-14 215.635484
1901-01-15 211.761290
1901-01-16 207.887097
1901-01-17 204.012903
1901-01-18 200.138710
1901-01-19 196.264516
1901-01-20 192.390323
1901-01-21 188.516129
1901-01-22 184.641935
1901-01-23 180.767742
1901-01-24 176.893548
1901-01-25 173.019355
1901-01-26 169.145161
1901-01-27 165.270968
1901-01-28 161.396774
1901-01-29 157.522581
1901-01-30 153.648387
1901-01-31 149.774194
1901-02-01 145.900000
Freq: D, Name: Sales, dtype: float64interpolated.plot()
pyplot.show()
pyplot.close()
## 7.4 Downsampling Data
# downsample to quarterly intervals
from pandas import read_csv
from datetime import datetime
from matplotlib import pyplot
def parser(x):
return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")<string>:1: FutureWarning: The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.resample = series.resample('Q')
quarterly_mean_sales = resample.mean()
print(quarterly_mean_sales.head())Month
1901-03-31 198.333333
1901-06-30 156.033333
1901-09-30 216.366667
1901-12-31 215.100000
1902-03-31 184.633333
Freq: Q-DEC, Name: Sales, dtype: float64quarterly_mean_sales.plot()
pyplot.show()
pyplot.close()
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")<string>:2: FutureWarning: The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.resample = series.resample('Y')
yearly_mean_sales = resample.sum()
print(yearly_mean_sales.head())Month
1901-12-31 2357.5
1902-12-31 3153.5
1903-12-31 5742.6
Freq: A-DEC, Name: Sales, dtype: float64yearly_mean_sales.plot()
pyplot.show()
pyplot.close()### Chapter 8 Power Transforms
## 8.1 Airline Passengers Dataset
# load and plot a time series
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(series)
# histogram
pyplot.subplot(212)
pyplot.hist(series)
pyplot.show()
pyplot.close()
## 8.2 Square Root Transform
# contrive a quadratic time series
from matplotlib import pyplot
series = [i**2 for i in range(1,100)]
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(series)
# histogram
pyplot.subplot(212)
pyplot.hist(series)
pyplot.show()
pyplot.close()
# square root transform a contrived quadratic time series
from matplotlib import pyplot
from numpy import sqrt
series = [i**2 for i in range(1,100)]
# sqrt transform
transform = series = sqrt(series)
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(transform)
# histogram
pyplot.subplot(212)
pyplot.hist(transform)
pyplot.show()
pyplot.close()
# square root transform a time series
from pandas import read_csv
from pandas import DataFrame
from numpy import sqrt
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
dataframe = DataFrame(series.values)
dataframe.columns = ['passengers']
dataframe['passengers'] = sqrt(dataframe['passengers'])
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(dataframe['passengers'])
# histogram
pyplot.subplot(212)
pyplot.hist(dataframe['passengers'])
pyplot.show()
pyplot.close()
## 8.3 Log Transform
# create and plot an exponential time series
from matplotlib import pyplot
from math import exp
series = [exp(i) for i in range(1,100)]
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(series)
# histogram
pyplot.subplot(212)
pyplot.hist(series)
pyplot.show()
pyplot.close()
# log transform a contrived exponential time series
from matplotlib import pyplot
from math import exp
from numpy import log
series = [exp(i) for i in range(1,100)]
transform = log(series)
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(transform)
# histogram
pyplot.subplot(212)
pyplot.hist(transform)
pyplot.show()
pyplot.close()
# log transform a time series
from pandas import read_csv
from pandas import DataFrame
from numpy import log
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
dataframe = DataFrame(series.values)
dataframe.columns = ['passengers']
dataframe['passengers'] = log(dataframe['passengers'])
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(dataframe['passengers'])
# histogram
pyplot.subplot(212)
pyplot.hist(dataframe['passengers'])
pyplot.show()
pyplot.close()
## 8.4 Box-Cox Transform
# manually box-cox transform a time series
from pandas import read_csv
from pandas import DataFrame
from scipy.stats import boxcox
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
dataframe = DataFrame(series.values)
dataframe.columns = ['passengers']
dataframe['passengers'] = boxcox(dataframe['passengers'], lmbda=0.0)
pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(dataframe['passengers'])
# histogram
pyplot.subplot(212)
pyplot.hist(dataframe['passengers'])
pyplot.show()
pyplot.close()
# automatically box-cox transform a time series
from pandas import read_csv
from pandas import DataFrame
from scipy.stats import boxcox
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
dataframe = DataFrame(series.values)
dataframe.columns = ['passengers']
dataframe['passengers'], lam = boxcox(dataframe['passengers'])
print('Lambda: %f' % lam)Lambda: 0.148023pyplot.figure(1)
# line plot
pyplot.subplot(211)
pyplot.plot(dataframe['passengers'])
# histogram
pyplot.subplot(212)
pyplot.hist(dataframe['passengers'])
pyplot.show()
pyplot.close()### Chapter 9 Moving Average Smoothing
## 9.4 Moving Average as Data Preparation
# moving average smoothing as data preparation
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# tail-rolling average transform
rolling = series.rolling(window=3)
rolling_mean = rolling.mean()
print(rolling_mean.head(10))Date
1959-01-01 NaN
1959-01-02 NaN
1959-01-03 32.333333
1959-01-04 31.000000
1959-01-05 35.000000
1959-01-06 34.666667
1959-01-07 39.333333
1959-01-08 39.000000
1959-01-09 42.000000
1959-01-10 36.000000
Name: Births, dtype: float64# plot original and transformed dataset
series.plot()
rolling_mean.plot(color='red')
pyplot.show()
pyplot.close()
# zoomed plot original and transformed dataset
series[:100].plot()
rolling_mean[:100].plot(color='red')
pyplot.show()
pyplot.close()
## 9.5 Moving Average as Feature Engineering
# moving average smoothing as feature engineering
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
df = DataFrame(series.values)
width = 3
lag1 = df.shift(1)
lag3 = df.shift(width - 1)
window = lag3.rolling(window=width)
means = window.mean()
dataframe = concat([means, lag1, df], axis=1)
dataframe.columns = ['mean', 't', 't+1']
print(dataframe.head(10)) mean t t+1
0 NaN NaN 35
1 NaN 35.0 32
2 NaN 32.0 30
3 NaN 30.0 31
4 32.333333 31.0 44
5 31.000000 44.0 29
6 35.000000 29.0 45
7 34.666667 45.0 43
8 39.333333 43.0 38
9 39.000000 38.0 27## 9.6 Moving Average as Prediction
# moving average smoothing as a forecast model
#!pip install scikit-learn
from math import sqrt
from pandas import read_csv
from numpy import mean
from sklearn.metrics import mean_squared_error
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# prepare situation
X = series.values
Xarray([35, 32, 30, 31, 44, 29, 45, 43, 38, 27, 38, 33, 55, 47, 45, 37, 50,
43, 41, 52, 34, 53, 39, 32, 37, 43, 39, 35, 44, 38, 24, 23, 31, 44,
38, 50, 38, 51, 31, 31, 51, 36, 45, 51, 34, 52, 47, 45, 46, 39, 48,
37, 35, 52, 42, 45, 39, 37, 30, 35, 28, 45, 34, 36, 50, 44, 39, 32,
39, 45, 43, 39, 31, 27, 30, 42, 46, 41, 36, 45, 46, 43, 38, 34, 35,
56, 36, 32, 50, 41, 39, 41, 47, 34, 36, 33, 35, 38, 38, 34, 53, 34,
34, 38, 35, 32, 42, 34, 46, 30, 46, 45, 54, 34, 37, 35, 40, 42, 58,
51, 32, 35, 38, 33, 39, 47, 38, 52, 30, 34, 40, 35, 42, 41, 42, 38,
24, 34, 43, 36, 55, 41, 45, 41, 37, 43, 39, 33, 43, 40, 38, 45, 46,
34, 35, 48, 51, 36, 33, 46, 42, 48, 34, 41, 35, 40, 34, 30, 36, 40,
39, 45, 38, 47, 33, 30, 42, 43, 41, 41, 59, 43, 45, 38, 37, 45, 42,
57, 46, 51, 41, 47, 26, 35, 44, 41, 42, 36, 45, 45, 45, 47, 38, 42,
35, 36, 39, 45, 43, 47, 36, 41, 50, 39, 41, 46, 64, 45, 34, 38, 44,
48, 46, 44, 37, 39, 44, 45, 33, 44, 38, 46, 46, 40, 39, 44, 48, 50,
41, 42, 51, 41, 44, 38, 68, 40, 42, 51, 44, 45, 36, 57, 44, 42, 53,
42, 34, 40, 56, 44, 53, 55, 39, 59, 55, 73, 55, 44, 43, 40, 47, 51,
56, 49, 54, 56, 47, 44, 43, 42, 45, 50, 48, 43, 40, 59, 41, 42, 51,
49, 45, 43, 42, 38, 47, 38, 36, 42, 35, 28, 44, 36, 45, 46, 48, 49,
43, 42, 59, 45, 52, 46, 42, 40, 40, 45, 35, 35, 40, 39, 33, 42, 47,
51, 44, 40, 57, 49, 45, 49, 51, 46, 44, 52, 45, 32, 46, 41, 34, 33,
36, 49, 43, 43, 34, 39, 35, 52, 47, 52, 39, 40, 42, 42, 53, 39, 40,
38, 44, 34, 37, 52, 48, 55, 50], dtype=int64)window = 3
history = [X[i] for i in range(window)]
history[35, 32, 30]test = [X[i] for i in range(window, len(X))]
test[31, 44, 29, 45, 43, 38, 27, 38, 33, 55, 47, 45, 37, 50, 43, 41, 52, 34, 53, 39, 32, 37, 43, 39, 35, 44, 38, 24, 23, 31, 44, 38, 50, 38, 51, 31, 31, 51, 36, 45, 51, 34, 52, 47, 45, 46, 39, 48, 37, 35, 52, 42, 45, 39, 37, 30, 35, 28, 45, 34, 36, 50, 44, 39, 32, 39, 45, 43, 39, 31, 27, 30, 42, 46, 41, 36, 45, 46, 43, 38, 34, 35, 56, 36, 32, 50, 41, 39, 41, 47, 34, 36, 33, 35, 38, 38, 34, 53, 34, 34, 38, 35, 32, 42, 34, 46, 30, 46, 45, 54, 34, 37, 35, 40, 42, 58, 51, 32, 35, 38, 33, 39, 47, 38, 52, 30, 34, 40, 35, 42, 41, 42, 38, 24, 34, 43, 36, 55, 41, 45, 41, 37, 43, 39, 33, 43, 40, 38, 45, 46, 34, 35, 48, 51, 36, 33, 46, 42, 48, 34, 41, 35, 40, 34, 30, 36, 40, 39, 45, 38, 47, 33, 30, 42, 43, 41, 41, 59, 43, 45, 38, 37, 45, 42, 57, 46, 51, 41, 47, 26, 35, 44, 41, 42, 36, 45, 45, 45, 47, 38, 42, 35, 36, 39, 45, 43, 47, 36, 41, 50, 39, 41, 46, 64, 45, 34, 38, 44, 48, 46, 44, 37, 39, 44, 45, 33, 44, 38, 46, 46, 40, 39, 44, 48, 50, 41, 42, 51, 41, 44, 38, 68, 40, 42, 51, 44, 45, 36, 57, 44, 42, 53, 42, 34, 40, 56, 44, 53, 55, 39, 59, 55, 73, 55, 44, 43, 40, 47, 51, 56, 49, 54, 56, 47, 44, 43, 42, 45, 50, 48, 43, 40, 59, 41, 42, 51, 49, 45, 43, 42, 38, 47, 38, 36, 42, 35, 28, 44, 36, 45, 46, 48, 49, 43, 42, 59, 45, 52, 46, 42, 40, 40, 45, 35, 35, 40, 39, 33, 42, 47, 51, 44, 40, 57, 49, 45, 49, 51, 46, 44, 52, 45, 32, 46, 41, 34, 33, 36, 49, 43, 43, 34, 39, 35, 52, 47, 52, 39, 40, 42, 42, 53, 39, 40, 38, 44, 34, 37, 52, 48, 55, 50]predictions = list()
# walk forward over time steps in test
for t in range(len(test)):
length = len(history)
yhat = mean([history[i] for i in range(length-window,length)])
obs = test[t]
predictions.append(yhat)
history.append(obs)
#print('predicted=%f, expected=%f' % (yhat, obs))
rmse = sqrt(mean_squared_error(test, predictions))
print('Test RMSE: %.3f' % rmse)Test RMSE: 7.834# plot
pyplot.plot(test)
pyplot.plot(predictions, color='red')
pyplot.show()
# zoom plot
pyplot.plot(test[:100])
pyplot.plot(predictions[:100], color='red')
pyplot.show()
### Chapter 10 A Gentle Introduction to White Noise
## 10.4 Example of White Noise Time Series
from random import gauss
from random import seed
from pandas import Series
from pandas.plotting import autocorrelation_plot
from matplotlib import pyplot
# seed random number generator
seed(1)
# create white noise series
series = [gauss(0.0, 1.0) for i in range(1000)]
series = Series(series)
# summary stats
print(series.describe())count 1000.000000
mean -0.013222
std 1.003685
min -2.961214
25% -0.684192
50% -0.010934
75% 0.703915
max 2.737260
dtype: float64# line plot
series.plot()
pyplot.show()
pyplot.close()
# histogram plot
series.hist()
pyplot.show()
pyplot.close()
# autocorrelation
autocorrelation_plot(series)
pyplot.show()
pyplot.close()
# calculate and plot a white noise series
from random import gauss
from random import seed
from pandas import Series
from pandas.plotting import autocorrelation_plot
from matplotlib import pyplot
# seed random number generator
seed(1)
# create white noise series
series = [gauss(0.0, 1.0) for i in range(1000)]
series = Series(series)
# summary stats
print(series.describe())count 1000.000000
mean -0.013222
std 1.003685
min -2.961214
25% -0.684192
50% -0.010934
75% 0.703915
max 2.737260
dtype: float64# line plot
series.plot()
pyplot.show()
pyplot.close()
# histogram plot
series.hist()
pyplot.show()
pyplot.close()
# autocorrelation
autocorrelation_plot(series)
pyplot.show()
pyplot.close()### Chapter 11 A Gentle Introduction to the Random Walk
#!pip install statsmodels
# 11.1 Random Series
# create and plot a random series
from random import seed
from random import randrange
from matplotlib import pyplot
seed(1)
series = [randrange(10) for i in range(1000)]
pyplot.plot(series)
pyplot.show()
pyplot.close()
# 11.2 Random Walk
# create and plot a random walk
from random import seed
from random import random
from matplotlib import pyplot
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1,1000):
movement = -1 if random() < 0.5 else 1
value = random_walk[i-1] + movement
random_walk.append(value)
pyplot.plot(random_walk)
pyplot.show()
pyplot.close()
# 11.3 Random Walk and Autocorrelation
# plot the autocorrelation of a random walk
from random import seed
from random import random
from matplotlib import pyplot
from pandas.plotting import autocorrelation_plot
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1000):
movement = -1 if random() < 0.5 else 1
value = random_walk[i-1] + movement
random_walk.append(value)
autocorrelation_plot(random_walk)
pyplot.show()
pyplot.close()
# 11.4 Random Walk and Stationarity
# calculate the stationarity of a random walk
from random import seed
from random import random
from statsmodels.tsa.stattools import adfuller
# generate random walk
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1,1000):
movement = -1 if random() < 0.5 else 1
value = random_walk[i-1] + movement
random_walk.append(value)
print()# statistical test
result = adfuller(random_walk)
print('ADF Statistic: %f' % result[0])ADF Statistic: 0.341605print('p-value: %f' % result[1])p-value: 0.979175print('Critical Values:')Critical Values:for key, value in result[4].items():
print('\t%s: %.3f' % (key, value)) 1%: -3.437
5%: -2.864
10%: -2.568# calculate and plot a differenced random walk
from random import seed
from random import random
from matplotlib import pyplot
# create random walk
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1000):
movement = -1 if random() < 0.5 else 1
value = random_walk[i-1] + movement
random_walk.append(value)
# take difference
diff = list()
for i in range(1, len(random_walk)):
value = random_walk[i] - random_walk[i - 1]
diff.append(value)
# line plot
pyplot.plot(diff)
pyplot.show()
pyplot.close()
# plot the autocorrelation of a differenced random walk
from random import seed
from random import random
from matplotlib import pyplot
from pandas.plotting import autocorrelation_plot
# create random walk
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1, 1000):
movement = -1 if random() < 0.5 else 1
value = random_walk[i-1] + movement
random_walk.append(value)
# take difference
diff = list()
for i in range(1, len(random_walk)):
value = random_walk[i] - random_walk[i - 1]
diff.append(value)
# line plot
autocorrelation_plot(diff)
pyplot.show()
pyplot.close()
# 11.5 Predicting a Random Walk
# persistence forecasts for a random walk
from random import seed
from random import random
from sklearn.metrics import mean_squared_error
from math import sqrt
# generate the random walk
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1, 1000):
movement = -1 if random() < 0.5 else 1
value = random_walk[i-1] + movement
random_walk.append(value)
# prepare dataset
train_size = int(len(random_walk) * 0.66)
train, test = random_walk[0:train_size], random_walk[train_size:]
# persistence
predictions = list()
history = train[-1]
for i in range(len(test)):
yhat = history
predictions.append(yhat)
history = test[i]
rmse = sqrt(mean_squared_error(test, predictions))
print('Persistence RMSE: %.3f' % rmse)Persistence RMSE: 1.000# random predictions for a random walk
from random import seed
from random import random
from sklearn.metrics import mean_squared_error
from math import sqrt
# generate the random walk
seed(1)
random_walk = list()
random_walk.append(-1 if random() < 0.5 else 1)
for i in range(1000):
movement = -1 if random() < 0.5 else 1
value = random_walk[i-1] + movement
random_walk.append(value)
# prepare dataset
train_size = int(len(random_walk) * 0.66)
train, test = random_walk[0:train_size], random_walk[train_size:]
# random prediction
predictions = list()
history = train[-1]
for i in range(len(test)):
yhat = history + (-1 if random() < 0.5 else 1)
predictions.append(yhat)
history = test[i]
rmse = sqrt(mean_squared_error(test, predictions))
print('Random RMSE: %.3f' % rmse)Random RMSE: 9.827### Chapter 12 Decompose Time Series Data
# 12.4 Automatic Time Series Decomposition
# 12.4.1 Additive Decomposition
# additive decompose a contrived additive time series
from random import randrange
from matplotlib import pyplot
from statsmodels.tsa.seasonal import seasonal_decompose
series = [i+randrange(10) for i in range(1,100)]
result = seasonal_decompose(series, model='additive', period=1)
result.plot()
pyplot.show()
pyplot.close()
# 12.4.2 Multiplicative Decomposition
# multiplicative decompose a contrived multiplicative time series
from matplotlib import pyplot
from statsmodels.tsa.seasonal import seasonal_decompose
series = [i**2.0 for i in range(1,100)]
result = seasonal_decompose(series, model='multiplicative', period=1)
result.plot()
pyplot.show()
pyplot.close()
# 12.4.3 Airline Passengers Dataset
# multiplicative decompose time series
from pandas import read_csv
from matplotlib import pyplot
from statsmodels.tsa.seasonal import seasonal_decompose
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
result = seasonal_decompose(series, model='multiplicative')
result.plot()
pyplot.show()
pyplot.close()### Chapter 13 Use and Remove Trends
# detrend a time series using differencing
from pandas import read_csv
from datetime import datetime
from matplotlib import pyplot
def parser(x):
return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")<string>:1: FutureWarning: The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.X = series.values
diff = list()
for i in range(1, len(X)):
value = X[i] - X[i - 1]
diff.append(value)
pyplot.plot(diff)
pyplot.show()
# use a linear model to detrend a time series
from pandas import read_csv
from datetime import datetime
from sklearn.linear_model import LinearRegression
from matplotlib import pyplot
import numpy
def parser(x):
return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")<string>:1: FutureWarning: The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.# fit linear model
X = [i for i in range(0, len(series))]
X = numpy.reshape(X, (len(X), 1))
y = series.values
model = LinearRegression()
model.fit(X, y)LinearRegression()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.
LinearRegression()
# calculate trend
trend = model.predict(X)
# plot trend
pyplot.plot(y)
pyplot.plot(trend)
pyplot.show()
# detrend
detrended = [y[i]-trend[i] for i in range(0, len(series))]
# plot detrended
pyplot.plot(detrended)
pyplot.show()
pyplot.close()### Chapter 14 Use and Remove Seasonality
# 14.3 Seasonal Adjustment with Differencing
# deseasonalize a time series using differencing
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
X = series.values
diff = list()
days_in_year = 365
for i in range(days_in_year, len(X)):
value = X[i] - X[i - days_in_year]
diff.append(value)
pyplot.plot(diff)
pyplot.show()
pyplot.close()
# calculate and plot monthly average
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
resample = series.resample('M')
monthly_mean = resample.mean()
print(monthly_mean.head(13))Date
1981-01-31 17.712903
1981-02-28 17.678571
1981-03-31 13.500000
1981-04-30 12.356667
1981-05-31 9.490323
1981-06-30 7.306667
1981-07-31 7.577419
1981-08-31 7.238710
1981-09-30 10.143333
1981-10-31 10.087097
1981-11-30 11.890000
1981-12-31 13.680645
1982-01-31 16.567742
Freq: M, Name: Temp, dtype: float64monthly_mean.plot()
pyplot.show()
pyplot.close()
# deseasonalize monthly data by differencing
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
resample = series.resample('M')
monthly_mean = resample.mean()
X = series.values
diff = list()
months_in_year = 12
for i in range(months_in_year, len(monthly_mean)):
value = monthly_mean[i] - monthly_mean[i - months_in_year]
diff.append(value)<string>:2: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`pyplot.plot(diff)
pyplot.show()
pyplot.close()
# deseasonalize a time series using month-based differencing
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
X = series.values
diff = list()
days_in_year = 365
for i in range(days_in_year, len(X)):
month_str = str(series.index[i].year-1)+'-'+str(series.index[i].month)
month_mean_last_year = series[month_str].mean()
value = X[i] - month_mean_last_year
diff.append(value)
pyplot.plot(diff)
pyplot.show()
pyplot.close()
# model seasonality with a polynomial model
from pandas import read_csv
from matplotlib import pyplot
from numpy import polyfit
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
# fit polynomial: x^2*b1 + x*b2 + ... + bn
X = [i%365 for i in range(0, len(series))]
y = series.values
degree = 4
coef = polyfit(X, y, degree)
print('Coefficients: %s' % coef)Coefficients: [-1.17308000e-08 9.30253946e-06 -2.15977594e-03 1.19147966e-01
1.38980178e+01]# create curve
curve = list()
for i in range(len(X)):
value = coef[-1]
for d in range(degree):
value += X[i]**(degree-d) * coef[d]
curve.append(value)
# plot curve over original data
pyplot.plot(series.values)
pyplot.plot(curve, color='red', linewidth=3)
pyplot.show()
pyplot.close()
# deseasonalize by differencing with a polynomial model
from pandas import read_csv
from matplotlib import pyplot
from numpy import polyfit
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
# fit polynomial: x^2*b1 + x*b2 + ... + bn
X = [i%365 for i in range(0, len(series))]
y = series.values
degree = 4
coef = polyfit(X, y, degree)
# create curve
curve = list()
for i in range(len(X)):
value = coef[-1]
for d in range(degree):
value += X[i]**(degree-d) * coef[d]
curve.append(value)
# create seasonally adjusted
values = series.values
diff = list()
for i in range(len(values)):
value = values[i] - curve[i]
diff.append(value)
pyplot.plot(diff)
pyplot.show()
pyplot.close()### Chapter 15 Stationarity in Time Series Data
## 15.1 Stationary Time Series
# load time series data
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
series.plot()
pyplot.show()
pyplot.close()
#15.2 Non-Stationary Time Series
# load time series data
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
series.plot()
pyplot.show()
pyplot.close()
# 15.6.1 Daily Births Dataset
# plot a histogram of a time series
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
series.hist()
pyplot.show()
pyplot.close()
# calculate statistics of partitioned time series data
from pandas import read_csv
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
split = int(len(X) / 2)
X1, X2 = X[0:split], X[split:]
mean1, mean2 = X1.mean(), X2.mean()
var1, var2 = X1.var(), X2.var()
print('mean1=%f, mean2=%f' % (mean1, mean2))mean1=39.763736, mean2=44.185792print('variance1=%f, variance2=%f' % (var1, var2))variance1=49.213410, variance2=48.708651# 15.6.2 Airline Passengers Dataset
# calculate statistics of partitioned time series data
from pandas import read_csv
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
split = int(len(X) / 2)
X1, X2 = X[0:split], X[split:]
mean1, mean2 = X1.mean(), X2.mean()
var1, var2 = X1.var(), X2.var()
print('mean1=%f, mean2=%f' % (mean1, mean2))mean1=182.902778, mean2=377.694444print('variance1=%f, variance2=%f' % (var1, var2))variance1=2244.087770, variance2=7367.962191# plot a histogram of a time series
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
series.hist()
pyplot.show()
# histogram and line plot of log transformed time series
from pandas import read_csv
from matplotlib import pyplot
from numpy import log
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
X = log(X)
pyplot.hist(X)
pyplot.show()
pyplot.plot(X)
pyplot.show()
# calculate statistics of partitioned log transformed time series data
from pandas import read_csv
from numpy import log
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
X = log(X)
split = int(len(X) / 2)
X1, X2 = X[0:split], X[split:]
mean1, mean2 = X1.mean(), X2.mean()
var1, var2 = X1.var(), X2.var()
print('mean1=%f, mean2=%f' % (mean1, mean2))mean1=5.175146, mean2=5.909206print('variance1=%f, variance2=%f' % (var1, var2))variance1=0.068375, variance2=0.049264# 15.7 Augmented Dickey-Fuller test
# calculate stationarity test of time series data
from pandas import read_csv
from statsmodels.tsa.stattools import adfuller
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
result = adfuller(X)
print('ADF Statistic: %f' % result[0])ADF Statistic: -4.808291print('p-value: %f' % result[1])p-value: 0.000052print('Critical Values:')Critical Values:for key, value in result[4].items():
print('\t%s: %.3f' % (key, value)) 1%: -3.449
5%: -2.870
10%: -2.571
# calculate stationarity test of time series data
from pandas import read_csv
from statsmodels.tsa.stattools import adfuller
series = read_csv('monthly-airline-passengers.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
result = adfuller(X)
print('ADF Statistic: %f' % result[0])ADF Statistic: 0.815369print('p-value: %f' % result[1])p-value: 0.991880print('Critical Values:')Critical Values:for key, value in result[4].items():
print('\t%s: %.3f' % (key, value)) 1%: -3.482
5%: -2.884
10%: -2.579### Chapter 16 Backtest Forecast Models
# calculate a train-test split of a time series dataset
from pandas import read_csv
series = read_csv('sunspots.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
train_size = int(len(X) * 0.66)
train, test = X[0:train_size], X[train_size:len(X)]
print('Observations: %d' % (len(X)))Observations: 2820print('Training Observations: %d' % (len(train)))Training Observations: 1861print('Testing Observations: %d' % (len(test)))Testing Observations: 959# plot train-test split of time series data
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('sunspots.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
train_size = int(len(X) * 0.66)
train, test = X[0:train_size], X[train_size:len(X)]
print('Observations: %d' % (len(X)))Observations: 2820print('Training Observations: %d' % (len(train)))Training Observations: 1861print('Testing Observations: %d' % (len(test)))Testing Observations: 959pyplot.plot(train)
pyplot.plot([None for i in train] + [x for x in test])
pyplot.show()
# calculate repeated train-test splits of time series data
from pandas import read_csv
from sklearn.model_selection import TimeSeriesSplit
from matplotlib import pyplot
series = read_csv('sunspots.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
splits = TimeSeriesSplit(n_splits=3)
pyplot.figure(1)
index = 1
for train_index, test_index in splits.split(X):
train = X[train_index]
test = X[test_index]
#print('Observations: %d' % (len(train) + len(test)))
#print('Training Observations: %d' % (len(train)))
#print('Testing Observations: %d' % (len(test)))
pyplot.subplot(310 + index)
pyplot.plot(train)
pyplot.plot([None for i in train] + [x for x in test])
index += 1
pyplot.show()
# 16.5 Walk Forward Validation
# walk forward evaluation model for time series data
from pandas import read_csv
series = read_csv('sunspots.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
X = series.values
n_train = 500
n_records = len(X)
for i in range(n_train, n_records):
train, test = X[0:i], X[i:i+1]
print('train=%d, test=%d' % (len(train), len(test)))train=2819, test=1### Chapter 17 Forecasting Performance Measures
# calculate forecast error
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
forecast_errors = [expected[i]-predictions[i] for i in range(len(expected))]
print('Forecast Errors: %s' % forecast_errors)Forecast Errors: [-0.2, 0.09999999999999998, -0.1, -0.09999999999999998, -0.2]# 17.2 Mean Forecast Error (or Forecast Bias)
# calculate mean forecast error
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
forecast_errors = [expected[i]-predictions[i] for i in range(len(expected))]
bias = sum(forecast_errors) * 1.0/len(expected)
print('Bias: %f' % bias)Bias: -0.100000# 17.3 Mean Absolute Error
# calculate mean absolute error
from sklearn.metrics import mean_absolute_error
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
mae = mean_absolute_error(expected, predictions)
print('MAE: %f' % mae)MAE: 0.140000# 17.4 Mean Squared Error
# calculate mean squared error
from sklearn.metrics import mean_squared_error
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
mse = mean_squared_error(expected, predictions)
print('MSE: %f' % mse)MSE: 0.022000#17.5 Root Mean Squared Error
# calculate root mean squared error
from sklearn.metrics import mean_squared_error
from math import sqrt
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
mse = mean_squared_error(expected, predictions)
rmse = sqrt(mse)
print('RMSE: %f' % rmse)RMSE: 0.148324# Chapter 18 Persistence Model for Forecasting
# 18.4.1 Step 1: Define the Supervised Learning Problem
# Create lagged dataset
# evaluate a persistence forecast model
from pandas import read_csv
from datetime import datetime
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
from sklearn.metrics import mean_squared_error
from math import sqrt
# load dataset
def parser(x):
return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('monthly-shampoo-sales.csv', header=0, index_col=0, parse_dates=True,date_parser=parser).squeeze("columns")<string>:1: FutureWarning: The argument 'date_parser' is deprecated and will be removed in a future version. Please use 'date_format' instead, or read your data in as 'object' dtype and then call 'to_datetime'.# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
print(dataframe.head(5)) t t+1
0 NaN 266.0
1 266.0 145.9
2 145.9 183.1
3 183.1 119.3
4 119.3 180.3# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
def model_persistence(x):
return x
# walk-forward validation
predictions = list()
for x in test_X:
yhat = model_persistence(x)
predictions.append(yhat)
rmse = sqrt(mean_squared_error(test_y, predictions))
print('Test RMSE: %.3f' % rmse)Test RMSE: 133.156# plot predictions and expected results
pyplot.plot(train_y)
pyplot.plot([None for i in train_y] + [x for x in test_y])
pyplot.plot([None for i in train_y] + [x for x in predictions])
pyplot.show()
### Chapter 19 Visualize Residual Forecast Errors
# create lagged dataset
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
# calculate residuals from a persistence forecast
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
residuals = DataFrame(residuals)
print(residuals.head()) 0
0 9.0
1 -10.0
2 3.0
3 -6.0
4 30.0# 19.4 Residual Line Plot
# line plot of residual errors
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
residuals = DataFrame(residuals)
# plot residuals
residuals.plot()
pyplot.show()
#19.5 Residual Summary Statistics
# summary statistics of residual errors
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
residuals = DataFrame(residuals)
# summary statistics
print(residuals.describe()) 0
count 125.000000
mean 0.064000
std 9.187776
min -28.000000
25% -6.000000
50% -1.000000
75% 5.000000
max 30.000000#19.6 Residual Histogram and Density Plots
# density plots of residual errors
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
residuals = DataFrame(residuals)
# histogram plot
residuals.hist()array([[<Axes: title={'center': '0'}>]], dtype=object)pyplot.show()
# density plot
residuals.plot(kind='kde')
pyplot.show()
# 19.7 Residual Q-Q Plot
# qq plot of residual errors
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
import numpy
from statsmodels.graphics.gofplots import qqplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
residuals = numpy.array(residuals)
qqplot(residuals, line='r')
pyplot.show()
# 19.8 Residual Autocorrelation Plot
# autoregression plot of residual errors
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
from pandas.plotting import autocorrelation_plot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0, parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
predictions = [x for x in test_X]
# calculate residuals
residuals = [test_y[i]-predictions[i] for i in range(len(predictions))]
residuals = DataFrame(residuals)
autocorrelation_plot(residuals)
pyplot.show()
### Chapter 20 Reframe Time Series Forecasting Problems
# reframe precision of regression forecast
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
# load data
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# round forecast to nearest 5
for i in range(len(dataframe['t+1'])):
dataframe.loc[i,'t+1'] = int(dataframe.loc[i,'t+1'] / 5) * 5.0
print(dataframe.head(5)) t t+1
0 NaN 20.0
1 20.7 15.0
2 17.9 15.0
3 18.8 10.0
4 14.6 15.0# 20.5 Classification Framings
# reframe regression as classification
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
# load data
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
# Create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# make discrete
for i in range(len(dataframe['t+1'])):
value = dataframe['t+1'][i]
if value < 10.0:
dataframe.loc[i,'t+1'] = 0
elif value >= 25.0:
dataframe.loc[i,'t+1'] = 2
else:
dataframe.loc[i,'t+1'] = 1
print(dataframe.head(5)) t t+1
0 NaN 1.0
1 20.7 1.0
2 17.9 1.0
3 18.8 1.0
4 14.6 1.0# 20.6 Time Horizon Framings
# reframe time horizon of forecast
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
# load data
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values, values.shift(-1),
values.shift(-2), values.shift(-3), values.shift(-4), values.shift(-5),
values.shift(-6)], axis=1)
dataframe.columns = ['t', 't+1', 't+2', 't+3', 't+4', 't+5', 't+6', 't+7']
print(dataframe.head(14)) t t+1 t+2 t+3 t+4 t+5 t+6 t+7
0 NaN 20.7 17.9 18.8 14.6 15.8 15.8 15.8
1 20.7 17.9 18.8 14.6 15.8 15.8 15.8 17.4
2 17.9 18.8 14.6 15.8 15.8 15.8 17.4 21.8
3 18.8 14.6 15.8 15.8 15.8 17.4 21.8 20.0
4 14.6 15.8 15.8 15.8 17.4 21.8 20.0 16.2
5 15.8 15.8 15.8 17.4 21.8 20.0 16.2 13.3
6 15.8 15.8 17.4 21.8 20.0 16.2 13.3 16.7
7 15.8 17.4 21.8 20.0 16.2 13.3 16.7 21.5
8 17.4 21.8 20.0 16.2 13.3 16.7 21.5 25.0
9 21.8 20.0 16.2 13.3 16.7 21.5 25.0 20.7
10 20.0 16.2 13.3 16.7 21.5 25.0 20.7 20.6
11 16.2 13.3 16.7 21.5 25.0 20.7 20.6 24.8
12 13.3 16.7 21.5 25.0 20.7 20.6 24.8 17.7
13 16.7 21.5 25.0 20.7 20.6 24.8 17.7 15.5### Chapter 22 Autoregression Models for Forecasting
# lag plot of time series
from pandas import read_csv
from matplotlib import pyplot
from pandas.plotting import lag_plot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
lag_plot(series)
pyplot.show()
# correlation of lag=1
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
result = dataframe.corr()
print(result) t t+1
t 1.00000 0.77487
t+1 0.77487 1.00000#22.5 Autocorrelation Plots
# autocorrelation plot of time series
from pandas import read_csv
from matplotlib import pyplot
from pandas.plotting import autocorrelation_plot
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
autocorrelation_plot(series)
pyplot.show()
# autocorrelation plot of time series
from pandas import read_csv
from matplotlib import pyplot
from statsmodels.graphics.tsaplots import plot_acf
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
plot_acf(series, lags=31)
pyplot.show()
# 22.6 Persistence Model
# evaluate a persistence model
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
from sklearn.metrics import mean_squared_error
from math import sqrt
# load dataset
series = read_csv('daily-min-temperatures (1).csv', header=0, index_col=0,
parse_dates=True).squeeze("columns")
# create lagged dataset
values = DataFrame(series.values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t', 't+1']
# split into train and test sets
X = dataframe.values
train, test = X[1:len(X)-7], X[len(X)-7:]
train_X, train_y = train[:,0], train[:,1]
test_X, test_y = test[:,0], test[:,1]
# persistence model
def model_persistence(x):
return x
# walk-forward validation
predictions = list()
for x in test_X:
yhat = model_persistence(x)
predictions.append(yhat)
rmse = sqrt(mean_squared_error(test_y, predictions))
print('Test RMSE: %.3f' % rmse)Test RMSE: 1.850# plot predictions vs expected
pyplot.plot(test_y)
pyplot.plot(predictions, color='red')
pyplot.show()
Time Series Forecasting Using Generative AI: Leveraging AI for Precision Forecasting
Generative AI is a subset of artificial intelligence because it utilizes AI techniques, such as machine learning and pattern recognition, to generate new content, like images and text; just as how a painter uses brushes to create art, GenAI uses algorithms to create new content, making it a specialized tool within the broader scope of AI. For example, ChatGPT, a GenAI tool, uses AI algorithms to generate human-like text responses, making it a subset of AI.
GenAI involves leveraging AI to generate novel content, such as text, images, music, audio, and videos, by employing machine learning algorithms to identify patterns and relationships within human-created content. These learned patterns are then used to create new content, effectively mimicking human creativity. The emergence of GenAI has significant implications for language teaching and learning, which plays a vital role in today’s globalized world. Language proficiency enables individuals to communicate effectively, express ideas clearly, and navigate diverse cultural contexts
Multilayer Perceptron (MLP)
An MLP is a neural network that has at least three layers: an input layer, a hidden layer, and an output layer. Each layer performs operations on the outputs of its preceding layer
#!pip install IPython
#!pip install pandas
# 2.2.1 Multilayer Perceptron in Action
import pandas as pd
import numpy as np
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from IPython.display import display, Markdown
import matplotlib.pyplot as plt
from neuralforecast import NeuralForecast
from neuralforecast.losses.pytorch import MQLoss,DistributionLoss
from neuralforecast.models import MLP
def calculate_error_metrics(actual, predicted, num_predictors=1):
# convert inputs are numpy arrays
actual = np.array(actual)
predicted = np.array(predicted)
# Number of observations
n = len(actual)
# Calculate MSE
mse = mean_squared_error(actual, predicted)
# Calculate RMSE
rmse = np.sqrt(mse)
# Calculate MAPE
mape = mean_absolute_percentage_error(actual, predicted)
# Calculate R-squared
r2 = r2_score(actual, predicted)
# Calculate Adjusted R-squared
adjusted_r2 = 1 - ((1 - r2) * (n - 1) / (n - num_predictors - 1))
print(f'MSE : {mse}')
print(f'RMSE : {rmse}')
print(f'MAPE : {mape}')
print(f'r2 : {r2}')
print(f'adjusted_r2 : {adjusted_r2}')
from neuralforecast.utils import AirPassengersDF
Y_df = AirPassengersDF
Y_df = Y_df.reset_index(drop=True)
Y_df.head() unique_id ds y
0 1.0 1949-01-31 112.0
1 1.0 1949-02-28 118.0
2 1.0 1949-03-31 132.0
3 1.0 1949-04-30 129.0
4 1.0 1949-05-31 121.0train_data = Y_df.head(132)
test_data = Y_df.tail(12)
horizon = 12
model = MLP(h=horizon, input_size=12,
loss=DistributionLoss(distribution='Normal',
level=[80, 90]),
scaler_type='robust',
learning_rate=1e-3,
max_steps=200,
val_check_steps=10,
early_stop_patience_steps=2)
#fcst = NeuralForecast(models=[model],freq='M')
#fcst.fit(df=train_data, val_size=12)
#Y_hat_df = fcst.predict()
#Y_hat_df.head()
calculate_error_metrics(test_data[['y']],Y_hat_df['MLP'])MSE : 547.1504516601562
RMSE : 23.39124733014801
MAPE : 0.03406834229826927
r2 : 0.9012269377708435
adjusted_r2 : 0.8913496315479279
train_data.set_index('ds',inplace =True)
Y_hat_df.set_index('ds',inplace =True)
test_data.set_index('ds',inplace =True)
plt.figure(figsize=(25, 15))
y_past = train_data["y"]
y_pred = Y_hat_df['MLP']
y_test = test_data["y"]
plt.plot(y_past, label="Past time series values")
plt.plot(y_pred, label="Forecast")
plt.plot(y_test, label="Actual time series values")
plt.title('AirPassengers Forecast', fontsize=10)
plt.ylabel('Monthly Passengers', fontsize=10)
plt.xlabel('Timestamp [t]', fontsize=10)
plt.legend();
Temporal Convolutional Networks (TCN)
import pandas as pd
import numpy as np
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from IPython.display import display, Markdown
import matplotlib.pyplot as plt
from neuralforecast import NeuralForecast
from neuralforecast.losses.pytorch import GMM, MQLoss, DistributionLoss
from neuralforecast.auto import TCN
from neuralforecast.tsdataset import TimeSeriesDataset
from ray import tune
from neuralforecast.utils import AirPassengersDF as Y_df
def calculate_error_metrics(actual, predicted, num_predictors=1):
# convert inputs are numpy arrays
actual = np.array(actual)
predicted = np.array(predicted)
# Number of observations
n = len(actual)
# Calculate MSE
mse = mean_squared_error(actual, predicted)
# Calculate RMSE
rmse = np.sqrt(mse)
# Calculate MAPE
mape = mean_absolute_percentage_error(actual, predicted)
# Calculate R-squared
r2 = r2_score(actual, predicted)
# Calculate Adjusted R-squared
adjusted_r2 = 1 - ((1 - r2) * (n - 1) / (n - num_predictors - 1))
print(f'MSE : {mse}')
print(f'RMSE : {rmse}')
print(f'MAPE : {mape}')
print(f'r2 : {r2}')
print(f'adjusted_r2 : {adjusted_r2}')
Y_train_df = Y_df[Y_df.ds<='1959-12-31']
Y_test_df = Y_df[Y_df.ds>'1959-12-31']
dataset, *_ = TimeSeriesDataset.from_df(Y_train_df)calculate_error_metrics(Y_test_df[['y']],y_hat[['TCN']])MSE : 396.6375732421875
RMSE : 19.915761929742672
MAPE : 0.03342359513044357
r2 : 0.9283979535102844
adjusted_r2 : 0.9212377488613128Y_train_df.set_index('ds',inplace =True)
Y_test_df.set_index('ds',inplace =True)
plt.figure(figsize=(25, 15))<Figure size 2500x1500 with 0 Axes>y_past = Y_train_df["y"]
y_pred = y_hat[['TCN']]
y_test = Y_test_df["y"]
plt.plot(y_past, label="Past time series values")[<matplotlib.lines.Line2D object at 0x0000022F948CB1F0>]plt.plot(y_pred, label="Forecast")[<matplotlib.lines.Line2D object at 0x0000022F948DE1D0>]plt.plot(y_test, label="Actual time series values")[<matplotlib.lines.Line2D object at 0x0000022F948DF250>]plt.title('AirPassengers Forecast')Text(0.5, 1.0, 'AirPassengers Forecast')plt.ylabel('Monthly Passengers')Text(0, 0.5, 'Monthly Passengers')plt.xlabel('Timestamp [t]')Text(0.5, 0, 'Timestamp [t]')plt.legend();
Bidirectional Temporal Convolutional Network (BiTCN)
Bidirectional temporal convolutional network (BiTCN) architecture is developed by using two temporal convolutional networks for forecasting. The first network, called the forward network, encodes future covariates of the time series. The second network, called the backward network, encodes past observations and covariates. This technique helps in preserving temporal information of sequence data. The parameters of the output distribution are jointly estimated using these two networks. It is computationally more efficient than RNN methods like LSTM.
import pandas as pd
import numpy as np
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from IPython.display import display, Markdown
import matplotlib.pyplot as plt
from neuralforecast import NeuralForecast
from neuralforecast.models import BiTCN
from ray import tune
from neuralforecast.losses.pytorch import GMM, DistributionLoss
from neuralforecast.tsdataset import TimeSeriesDataset
from neuralforecast.utils import AirPassengersDF as Y_df
def calculate_error_metrics(actual, predicted, num_predictors=1):
# convert inputs are numpy arrays
actual = np.array(actual)
predicted = np.array(predicted)
# Number of observations
n = len(actual)
# Calculate MSE
mse = mean_squared_error(actual, predicted)
# Calculate RMSE
rmse = np.sqrt(mse)
# Calculate MAPE
mape = mean_absolute_percentage_error(actual, predicted)
# Calculate R-squared
r2 = r2_score(actual, predicted)
# Calculate Adjusted R-squared
adjusted_r2 = 1 - ((1 - r2) * (n - 1) / (n - num_predictors - 1))
print(f'MSE : {mse}')
print(f'RMSE : {rmse}')
print(f'MAPE : {mape}')
print(f'r2 : {r2}')
print(f'adjusted_r2 : {adjusted_r2}')
Y_train_df = Y_df[Y_df.ds<='1959-12-31']
Y_test_df = Y_df[Y_df.ds>'1959-12-31']
dataset, *_ = TimeSeriesDataset.from_df(Y_train_df)y_hat.set_index('ds',inplace =True)
y_hat.head() unique_id BiTCN ... BiTCN-mu-7 BiTCN-std-7
ds ...
1960-01-31 1.0 439.786194 ... 481.688324 57.548355
1960-02-29 1.0 458.564697 ... 485.977234 61.890015
1960-03-31 1.0 476.137177 ... 542.925232 53.607937
1960-04-30 1.0 485.645538 ... 522.366516 61.040913
1960-05-31 1.0 487.014954 ... 533.534668 54.934593
[5 rows x 21 columns]calculate_error_metrics(Y_test_df[['y']],y_hat[['BiTCN']])MSE : 3626.825439453125
RMSE : 60.2231304355156
MAPE : 0.10271913558244705
r2 : 0.3452759385108948
adjusted_r2 : 0.2798035323619843Y_train_df.set_index('ds',inplace =True)
Y_test_df.set_index('ds',inplace =True)
plt.figure(figsize=(25, 15))<Figure size 2500x1500 with 0 Axes>y_past = Y_train_df["y"]
y_pred = y_hat[['BiTCN']]
y_test = Y_test_df["y"]
plt.plot(y_past, label="Past time series values")[<matplotlib.lines.Line2D object at 0x0000022F9490A290>]plt.plot(y_pred, label="Forecast")[<matplotlib.lines.Line2D object at 0x0000022F9492CFA0>]plt.plot(y_test, label="Actual time series values")[<matplotlib.lines.Line2D object at 0x0000022F9492E020>]plt.title('AirPassengers Forecast', fontsize=30)Text(0.5, 1.0, 'AirPassengers Forecast')plt.ylabel('Monthly Passengers', fontsize=30)Text(0, 0.5, 'Monthly Passengers')plt.xlabel('Timestamp [t]', fontsize=30)Text(0.5, 0, 'Timestamp [t]')plt.legend();
Recurrent Neural Network (RNN)
In the RNN architecture, the output from the previous step becomes input for the current step. So, no points for guessing why sequence is important for RNN architecture. This is in contrast to the traditional neural networks where the inputs and outputs from each layer are independent of each other. The emergence of RNN architecture was majorly due to NLP (natural language processing) use cases that required prediction of the next word. RNN provided a solution to this problem by making use of a hidden layer. This layer works as memory, as it is used to remember some information about the sequence.
import pandas as pd
import numpy as np
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from IPython.display import display, Markdown
import matplotlib.pyplot as plt
from neuralforecast import NeuralForecast
from neuralforecast.losses.pytorch import GMM, MQLoss,DistributionLoss
from neuralforecast.models import RNN
from neuralforecast.tsdataset import TimeSeriesDataset
from ray import tune
from neuralforecast.utils import AirPassengersDF as Y_df
def calculate_error_metrics(actual, predicted, num_predictors=1):
# convert inputs are numpy arrays
actual = np.array(actual)
predicted = np.array(predicted)
# Number of observations
n = len(actual)
# Calculate MSE
mse = mean_squared_error(actual, predicted)
# Calculate RMSE
rmse = np.sqrt(mse)
# Calculate MAPE
mape = mean_absolute_percentage_error(actual, predicted)
# Calculate R-squared
r2 = r2_score(actual, predicted)
# Calculate Adjusted R-squared
adjusted_r2 = 1 - ((1 - r2) * (n - 1) / (n - num_predictors - 1))
print(f'MSE : {mse}')
print(f'RMSE : {rmse}')
print(f'MAPE : {mape}')
print(f'r2 : {r2}')
print(f'adjusted_r2 : {adjusted_r2}')
Y_train_df = Y_df[Y_df.ds<='1959-12-31']
Y_test_df = Y_df[Y_df.ds>'1959-12-31']
dataset, *_ = TimeSeriesDataset.from_df(Y_train_df)y_hat.set_index('ds',inplace =True)
y_hat.head() unique_id RNN-median ... RNN-hi-80 RNN-hi-90
ds ...
1960-01-31 1.0 424.830475 ... 454.812195 465.117096
1960-02-29 1.0 404.263458 ... 430.056854 443.181305
1960-03-31 1.0 461.973358 ... 494.343323 505.895691
1960-04-30 1.0 460.056946 ... 489.093109 502.045746
1960-05-31 1.0 477.102234 ... 507.555481 520.431641
[5 rows x 6 columns]calculate_error_metrics(Y_test_df[['y']],y_hat[['RNN-median']])MSE : 229.1974334716797
RMSE : 15.139267930507065
MAPE : 0.02436511404812336
r2 : 0.9586246609687805
adjusted_r2 : 0.9544871270656585Y_train_df.set_index('ds',inplace =True)
Y_test_df.set_index('ds',inplace =True)
plt.figure(figsize=(25, 15))<Figure size 2500x1500 with 0 Axes>y_past = Y_train_df["y"]
y_pred = y_hat[['RNN-median']]
y_test = Y_test_df["y"]
plt.plot(y_past)[<matplotlib.lines.Line2D object at 0x0000022F94A536A0>]plt.plot(y_pred)[<matplotlib.lines.Line2D object at 0x0000022F94A940A0>]plt.plot(y_test)[<matplotlib.lines.Line2D object at 0x0000022F94A95630>]
Long Short-Term Memory (LSTM)
The long short-term memory (LSTM) architecture is a modification to the RNN architecture to allow additional signal paths. These additional paths help in bypassing many processing steps encountered at each stage of the network. This modification helps in remembering information over a large number of time steps. While this modification improves the performance compared to the RNN, it also introduces additional complexity. This complexity has an impact on training speed compared to RNNs. Popular tools like Apple’s Siri and Google’s AlphaGo were based on LSTM.
import pandas as pd
import numpy as np
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from IPython.display import display, Markdown
import matplotlib.pyplot as plt
from neuralforecast import NeuralForecast
from neuralforecast.models import LSTM
from neuralforecast.tsdataset import TimeSeriesDataset
from neuralforecast.losses.pytorch import GMM, MQLoss,DistributionLoss
from neuralforecast.utils import AirPassengersDF as Y_df
from ray import tune
from neuralforecast.utils import AirPassengersDF as Y_df
Y_train_df = Y_df[Y_df.ds<='1959-12-31'] # 132 train data
Y_test_df = Y_df[Y_df.ds>'1959-12-31'] # 12 test data
dataset, *_ = TimeSeriesDataset.from_df(Y_train_df)#model.fit(dataset=dataset)
y_hat = fcst.predict()
Predicting: | | 0/? [00:00<?, ?it/s]
Predicting: 0%| | 0/1 [00:00<?, ?it/s]
Predicting DataLoader 0: 0%| | 0/1 [00:00<?, ?it/s]
Predicting DataLoader 0: 100%|##########| 1/1 [00:00<00:00, 101.07it/s]
Predicting DataLoader 0: 100%|##########| 1/1 [00:00<00:00, 101.07it/s]y_hat.set_index('ds',inplace =True)
y_hat.head() unique_id LSTM ... LSTM-hi-80 LSTM-hi-90
ds ...
1960-01-31 1.0 383.367157 ... 415.843842 425.110443
1960-02-29 1.0 384.142822 ... 423.718140 436.733459
1960-03-31 1.0 411.430298 ... 454.735565 468.883850
1960-04-30 1.0 432.984985 ... 476.166199 487.669250
1960-05-31 1.0 467.594727 ... 507.642029 516.355347
[5 rows x 7 columns]calculate_error_metrics(Y_test_df[['y']],y_hat[['LSTM']])MSE : 282.4617919921875
RMSE : 16.806599655855063
MAPE : 0.028130965307354927
r2 : 0.9490092396736145
adjusted_r2 : 0.9439101636409759Y_train_df.set_index('ds',inplace =True)
Y_test_df.set_index('ds',inplace =True)
plt.figure(figsize=(25, 15))<Figure size 2500x1500 with 0 Axes>y_past = Y_train_df["y"]
y_pred = y_hat[['LSTM']]
y_test = Y_test_df["y"]
plt.plot(y_past, label="Past time series values")[<matplotlib.lines.Line2D object at 0x0000022F94909BD0>]plt.plot(y_pred, label="Forecast")[<matplotlib.lines.Line2D object at 0x0000022F94AA6500>]plt.plot(y_test, label="Actual time series values")[<matplotlib.lines.Line2D object at 0x0000022F94AA6C20>]plt.title('AirPassengers Forecast', fontsize=30)Text(0.5, 1.0, 'AirPassengers Forecast')plt.ylabel('Monthly Passengers', fontsize=30)Text(0, 0.5, 'Monthly Passengers')plt.xlabel('Timestamp [t]', fontsize=30)Text(0.5, 0, 'Timestamp [t]')plt.legend();
Neural Networks Based on Autoregression (DeepAR)
The DeepAR model has some key benefits compared to traditional approaches. The major advantages that set DeepAR apart are as follows: (a) Relatively much less effort and time need to be spent on feature engineering to capture complex and group-dependent behavior. This is because the model learns seasonality and dependencies on given covariates across the time series. (b) It has an ability to provide forecasts for products with little to no historical data. This is because of learning from historical data of similar events. The DeepAR model has properties that help produce better forecasts by learning from historical behavior of all the time series taken together. (a) It incorporates a wide range of likelihood functions. This allows the time series modeling team to choose suitable functions based on statistical properties of the data. (b) The probabilistic forecasts are generated in the form of Monte Carlo Samples. These can be used to compute quantile estimates belonging to subranges in the prediction horizon. This is important as discussed earlier in this section, where we pointed out advantages of forecasts with probabilities compared to a point forecast value.
import pandas as pd
import numpy as np
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from IPython.display import display, Markdown
import matplotlib.pyplot as plt
from neuralforecast import NeuralForecast
from neuralforecast.losses.pytorch import MQLoss,DistributionLoss, GMM, PMM
from neuralforecast.tsdataset import TimeSeriesDataset
import pandas as pd
import pytorch_lightning as pl
import matplotlib.pyplot as plt
from neuralforecast.models import DeepAR
from neuralforecast.losses.pytorch import DistributionLoss, HuberMQLoss
from neuralforecast.utils import AirPassengers, AirPassengersPanel, AirPassengersStatic
from neuralforecast.utils import AirPassengers, AirPassengersPanel, AirPassengersStatic, AirPassengersPanel, AirPassengersStatic
AirPassengersPanel.head() unique_id ds y trend y_[lag12]
0 Airline1 1949-01-31 112.0 0 112.0
1 Airline1 1949-02-28 118.0 1 118.0
2 Airline1 1949-03-31 132.0 2 132.0
3 Airline1 1949-04-30 129.0 3 129.0
4 Airline1 1949-05-31 121.0 4 121.0print(AirPassengersStatic) unique_id airline1 airline2
0 Airline1 0 1
1 Airline2 1 0Y_train_df = AirPassengersPanel[AirPassengersPanel.ds<AirPassengersPanel['ds'].values[-12]]
Y_test_df = AirPassengersPanel[AirPassengersPanel.ds>=AirPassengersPanel['ds'].values[-12]].reset_index(drop=True)calculate_error_metrics(Y_test_df[['y']],Y_hat_df[['DeepAR']])MSE : 601.0615844726562
RMSE : 24.51655735360608
MAPE : 0.03668348118662834
r2 : 0.9785637259483337
adjusted_r2 : 0.9775893498550762Y_hat_df =Y_hat_df.reset_index(drop=False).drop(columns=['unique_id','ds'])
plot_df = pd.concat([Y_test_df, Y_hat_df], axis=1)
plot_df = pd.concat([Y_train_df, plot_df])
plt.figure(figsize=(20, 15))<Figure size 2000x1500 with 0 Axes>plot_df = plot_df[plot_df.unique_id=='Airline1'].drop('unique_id', axis=1)
plt.plot(plot_df['ds'], plot_df['y'], c='black', label='True')[<matplotlib.lines.Line2D object at 0x0000022F94A5E890>]plt.plot(plot_df['ds'], plot_df['DeepAR-median'], c='blue',label='median')[<matplotlib.lines.Line2D object at 0x0000022F94ADB400>]plt.fill_between(x=plot_df['ds'][-12:],
y1=plot_df['DeepAR-lo-90'][-12:].values,
y2=plot_df['DeepAR-hi-90'][-12:].values,
alpha=0.4, label='level 90')<matplotlib.collections.FillBetweenPolyCollection object at 0x0000022F94AD8D30>plt.title('AirPassengers Forecast', fontsize=10)Text(0.5, 1.0, 'AirPassengers Forecast')plt.ylabel('Monthly Passengers', fontsize=10)Text(0, 0.5, 'Monthly Passengers')plt.xlabel('Timestamp [t]', fontsize=10)Text(0.5, 0, 'Timestamp [t]')plt.legend()<matplotlib.legend.Legend object at 0x0000022F94A97160>plt.grid()
plt.plot()[]
Neural Basis Expansion Analysis (NBEATS)
NBEATS architecture was developed in the process of exploring the use of deep learning to solve univariate time series forecasting use cases. This architecture is designed using a deep learning network with multiple fully connected layers and relies on a network of backward and forward residual links. The intention of including this model in our journey to understand time series forecasting with GenAI is that this model was the first pure deep learning approach that performed better than existing statistical approaches in the Makridakis M-competition. The NBEATS model surpassed the winning solution of the M4 competition.
import pandas as pd
import numpy as np
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from IPython.display import display, Markdown
import matplotlib.pyplot as plt
from ray import tune
from neuralforecast.core import NeuralForecast
from neuralforecast.models import NBEATS, NHITS
from neuralforecast.utils import AirPassengersDF
Y_df = AirPassengersDF
Y_df = Y_df.reset_index(drop=True)
Y_df.head() unique_id ds y
0 1.0 1949-01-31 112.0
1 1.0 1949-02-28 118.0
2 1.0 1949-03-31 132.0
3 1.0 1949-04-30 129.0
4 1.0 1949-05-31 121.0train_data = Y_df.head(132)
test_data = Y_df.tail(12)Y_hat_df.head() index unique_id ds NBEATS
0 0 1.0 1960-01-31 421.764130
1 1 1.0 1960-02-29 437.824432
2 2 1.0 1960-03-31 454.093109
3 3 1.0 1960-04-30 454.956543
4 4 1.0 1960-05-31 502.554260calculate_error_metrics(test_data[['y']],Y_hat_df['NBEATS'])MSE : 984.4974975585938
RMSE : 31.37670310212011
MAPE : 0.06071820482611656
r2 : 0.8222759366035461
adjusted_r2 : 0.8045035302639008train_data.set_index('ds',inplace =True)
Y_hat_df.set_index('ds',inplace =True)
test_data.set_index('ds',inplace =True)
plt.figure(figsize=(25, 20))<Figure size 2500x2000 with 0 Axes>item_id = "airline_1"
y_past = train_data["y"]
y_pred = Y_hat_df['NBEATS']
y_test = test_data["y"]
plt.plot(y_past, label="Past time series values")[<matplotlib.lines.Line2D object at 0x0000022F94A945B0>]plt.plot(y_pred, label="Mean forecast")[<matplotlib.lines.Line2D object at 0x0000022F94C1F3A0>]plt.plot(y_test, label="Actual time series values")[<matplotlib.lines.Line2D object at 0x0000022F94C0C6D0>]plt.title('AirPassengers Forecast', fontsize=10)Text(0.5, 1.0, 'AirPassengers Forecast')plt.ylabel('Monthly Passengers', fontsize=10)Text(0, 0.5, 'Monthly Passengers')plt.xlabel('Timestamp [t]', fontsize=10)Text(0.5, 0, 'Timestamp [t]')plt.legend();
Transformers for Time Series
Transformers initially revolutionized natural language processing and have increasingly found their application in other realms such as computer vision, audio processing, bioinformatics, finance, robotics, and time series analysis.
Transformers are designed using stacked self-attention and pointwise, fully connected layers for both the encoder and decoder. Transformers are built on encoder-decoder architecture. The encoder applies the mathematical function to the data and transforms input to a certain representation, while the decoder applies the inverse function to recover back the original data.
Vanilla Transformer
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from neuralforecast import NeuralForecast
from neuralforecast.models import MLP,VanillaTransformer
from neuralforecast.losses.pytorch import MQLoss, DistributionLoss,MAE
from neuralforecast.tsdataset import TimeSeriesDataset
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from neuralforecast.utils import AirPassengers, AirPassengersPanel, AirPassengersStatic, augment_calendar_df
from neuralforecast.utils import AirPassengersDF as Y_df
print(Y_df) unique_id ds y
0 1.0 1949-01-31 112.0
1 1.0 1949-02-28 118.0
2 1.0 1949-03-31 132.0
3 1.0 1949-04-30 129.0
4 1.0 1949-05-31 121.0
.. ... ... ...
139 1.0 1960-08-31 606.0
140 1.0 1960-09-30 508.0
141 1.0 1960-10-31 461.0
142 1.0 1960-11-30 390.0
143 1.0 1960-12-31 432.0
[144 rows x 3 columns]Y_train_df = Y_df[Y_df.ds<='1959-12-31']
Y_test_df = Y_df[Y_df.ds>'1959-12-31']forecasts.head() unique_id ds VanillaTransformer
0 1.0 1960-01-31 390.566284
1 1.0 1960-02-29 410.820770
2 1.0 1960-03-31 431.629608
3 1.0 1960-04-30 442.306793
4 1.0 1960-05-31 468.161926calculate_error_metrics(Y_test_df[['y']],forecasts['VanillaTransformer'])MSE : 4713.275390625
RMSE : 68.65329846864607
MAPE : 0.10067308694124222
r2 : 0.1491471529006958
adjusted_r2 : 0.06406186819076543Y_train_df.set_index('ds',inplace =True)
forecasts.set_index('ds',inplace =True)
Y_test_df.set_index('ds',inplace =True)
plt.figure(figsize=(25, 20))<Figure size 2500x2000 with 0 Axes>y_past = Y_train_df["y"]
y_pred = forecasts['VanillaTransformer']
y_test = Y_test_df["y"]
plt.plot(y_past, label="Past time series values")[<matplotlib.lines.Line2D object at 0x0000022F94D88820>]plt.plot(forecasts, label="Forecast")[<matplotlib.lines.Line2D object at 0x0000022F94D8AE90>, <matplotlib.lines.Line2D object at 0x0000022F94D8AE60>]plt.plot(y_test, label="Actual time series values")[<matplotlib.lines.Line2D object at 0x0000022F94D78D30>]plt.title('AirPassengers Forecast', fontsize=10)Text(0.5, 1.0, 'AirPassengers Forecast')plt.ylabel('Monthly Passengers', fontsize=10)Text(0, 0.5, 'Monthly Passengers')plt.xlabel('Timestamp [t]', fontsize=10)Text(0.5, 0, 'Timestamp [t]')plt.legend();
Inverted Transformers
iTransformers take the transformer architecture but apply the attention and feedforward network on the inverted dimensions. This means that the time points of each individual series are embedded into variate tokens which can be used by the attention mechanisms to capture multivariate correlation, and the feedforward network learns nonlinear relationships.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from neuralforecast import NeuralForecast
from neuralforecast.models import iTransformer
from neuralforecast.losses.pytorch import MQLoss, DistributionLoss,MSE
from neuralforecast.tsdataset import TimeSeriesDataset
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from neuralforecast.utils import AirPassengersDF as Y_df
Y_train_df = Y_df[Y_df.ds<='1959-12-31']
Y_test_df = Y_df[Y_df.ds>'1959-12-31']
dataset, *_ = TimeSeriesDataset.from_df(Y_train_df)Y_test_df.head() unique_id ds y iTransformers
132 1.0 1960-01-31 417.0 431.639618
133 1.0 1960-02-29 391.0 395.834778
134 1.0 1960-03-31 419.0 450.014862
135 1.0 1960-04-30 461.0 444.223022
136 1.0 1960-05-31 472.0 475.579132calculate_error_metrics(Y_test_df[['y']],Y_test_df['iTransformers'])MSE : 494.6007995605469
RMSE : 22.23962228907107
MAPE : 0.035653989762067795
r2 : 0.9107133746147156
adjusted_r2 : 0.9017847120761872Y_train_df.set_index('ds',inplace =True)
Y_test_df.set_index('ds',inplace =True)
plt.figure(figsize=(25, 20))<Figure size 2500x2000 with 0 Axes>y_past = Y_train_df["y"]
y_pred = Y_test_df['iTransformers']
y_test = Y_test_df["y"]
plt.plot(y_past, label="Past time series values")[<matplotlib.lines.Line2D object at 0x0000022F94E21750>]plt.plot(y_pred, label="Forecast")[<matplotlib.lines.Line2D object at 0x0000022F94E9E740>]plt.plot(y_test, label="Actual time series values")[<matplotlib.lines.Line2D object at 0x0000022F94E9E470>]plt.title('AirPassengers Forecast', fontsize=10)Text(0.5, 1.0, 'AirPassengers Forecast')plt.ylabel('Monthly Passengers', fontsize=10)Text(0, 0.5, 'Monthly Passengers')plt.xlabel('Timestamp [t]', fontsize=10)Text(0.5, 0, 'Timestamp [t]')plt.legend();
DLinear
DLinear is a combination of a decomposition scheme used in Autoformer and FEDformer with linear layers. It decomposes raw data input into a trend component by a moving average kernel and a remainder or seasonal component. Then, two one-layer linear layers are applied to each component, and we sum up the two features to get the final prediction. By explicitly handling trends, DLinear enhances the performance of a vanilla linear when there is a clear trend in the data.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from neuralforecast import NeuralForecast
from neuralforecast.models import MLP,DLinear
from neuralforecast.losses.pytorch import MQLoss, DistributionLoss,MAE
from neuralforecast.tsdataset import TimeSeriesDataset
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from neuralforecast.utils import AirPassengers, AirPassengersPanel, AirPassengersStatic, augment_calendar_df
from neuralforecast.utils import AirPassengersDF
Y_df = AirPassengersDF
Y_df = Y_df.reset_index(drop=True)
Y_train_df = Y_df[Y_df.ds<='1959-12-31']
Y_test_df = Y_df[Y_df.ds>'1959-12-31']calculate_error_metrics(Y_test_df[['y']],forecasts['DLinear'])MSE : 1411.2696533203125
RMSE : 37.56686909126594
MAPE : 0.07227641344070435
r2 : 0.7452338933944702
adjusted_r2 : 0.7197572827339173Y_train_df.set_index('ds',inplace =True)
forecasts.set_index('ds',inplace =True)
Y_test_df.set_index('ds',inplace =True)
plt.figure(figsize=(25, 20))<Figure size 2500x2000 with 0 Axes>y_past = Y_train_df["y"]
y_pred = forecasts['DLinear']
y_test = Y_test_df["y"]
plt.plot(y_past, label="Past time series values")[<matplotlib.lines.Line2D object at 0x0000022F94B65E40>]plt.plot(y_pred, label="Forecast")[<matplotlib.lines.Line2D object at 0x0000022F94A8C460>]plt.plot(y_test, label="Actual time series values")[<matplotlib.lines.Line2D object at 0x0000022F94A8FAC0>]plt.title('AirPassengers Forecast', fontsize=10)Text(0.5, 1.0, 'AirPassengers Forecast')plt.ylabel('Monthly Passengers', fontsize=10)Text(0, 0.5, 'Monthly Passengers')plt.xlabel('Timestamp [t]', fontsize=10)Text(0.5, 0, 'Timestamp [t]')plt.legend();
NLinear
NLinear is part of the LTSF-Linear family of models specifically designed to boost the performance of Linear when there is a distribution shift in the dataset. NLinear first subtracts the input by the last value of the sequence, then the input goes through a linear layer, and the subtracted part is added back before making the final prediction. The subtraction and addition in NLinear are a simple normalization for the input sequence.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from neuralforecast import NeuralForecast
from neuralforecast.models import MLP,NLinear
from neuralforecast.losses.pytorch import MQLoss, DistributionLoss,MAE
from neuralforecast.tsdataset import TimeSeriesDataset
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from neuralforecast.utils import AirPassengers, AirPassengersPanel, AirPassengersStatic, augment_calendar_df
from neuralforecast.utils import AirPassengersDF
def calculate_error_metrics(actual, predicted, num_predictors=1):
# convert inputs are numpy arrays
actual = np.array(actual)
predicted = np.array(predicted)
# Number of observations
n = len(actual)
# Calculate MSE
mse = mean_squared_error(actual, predicted)
# Calculate RMSE
rmse = np.sqrt(mse)
# Calculate MAPE
mape = mean_absolute_percentage_error(actual, predicted)
# Calculate R-squared
r2 = r2_score(actual, predicted)
# Calculate Adjusted R-squared
adjusted_r2 = 1 - ((1 - r2) * (n - 1) / (n - num_predictors - 1))
print(f'MSE : {mse}')
print(f'RMSE : {rmse}')
print(f'MAPE : {mape}')
print(f'r2 : {r2}')
print(f'adjusted_r2 : {adjusted_r2}')
Y_df = AirPassengersDF
Y_df = Y_df.reset_index(drop=True)
Y_df.head() unique_id ds y
0 1.0 1949-01-31 112.0
1 1.0 1949-02-28 118.0
2 1.0 1949-03-31 132.0
3 1.0 1949-04-30 129.0
4 1.0 1949-05-31 121.0Y_train_df = Y_df[Y_df.ds<='1959-12-31']
Y_test_df = Y_df[Y_df.ds>'1959-12-31']calculate_error_metrics(Y_test_df[['y']],forecasts['NLinear'])MSE : 1214.2550048828125
RMSE : 34.846161982100874
MAPE : 0.06048479303717613
r2 : 0.7807995080947876
adjusted_r2 : 0.7588794589042663Y_train_df.set_index('ds',inplace =True)
forecasts.set_index('ds',inplace =True)
Y_test_df.set_index('ds',inplace =True)
plt.figure(figsize=(25, 20))<Figure size 2500x2000 with 0 Axes>y_past = Y_train_df["y"]
y_pred = forecasts['NLinear']
y_test = Y_test_df["y"]
plt.plot(y_past, label="Past time series values")[<matplotlib.lines.Line2D object at 0x0000022F94ADB820>]plt.plot(y_pred, label="Forecast")[<matplotlib.lines.Line2D object at 0x0000022F949D9300>]plt.plot(y_test, label="Actual time series values")[<matplotlib.lines.Line2D object at 0x0000022F94A8FB20>]plt.title('AirPassengers Forecast', fontsize=10)Text(0.5, 1.0, 'AirPassengers Forecast')plt.ylabel('Monthly Passengers', fontsize=10)Text(0, 0.5, 'Monthly Passengers')plt.xlabel('Timestamp [t]', fontsize=10)Text(0.5, 0, 'Timestamp [t]')plt.legend();
Patch Time Series Transformer (PatchTST)
PatchTST supports multivariate time series forecasting and self-supervised representation learning. It is based on the segmentation of time series into subseries-level patches, which serve as input tokens to the transformer. Channel independence is a property of PatchTST where each channel contains a single univariate time series. Channel independence helps share the same embedding and transformer weights across all the series. This helps the PatchTST model to apply attention weights separately to each channel, which helps in better capturing the unique features and patterns in each channel.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from neuralforecast import NeuralForecast
from neuralforecast.models import PatchTST
from neuralforecast.losses.pytorch import MQLoss, DistributionLoss,MAE
from neuralforecast.tsdataset import TimeSeriesDataset
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from neuralforecast.utils import AirPassengers, AirPassengersPanel, AirPassengersStatic, augment_calendar_df
from neuralforecast.utils import AirPassengersDF as Y_df
Y_train_df = Y_df[Y_df.ds<='1959-12-31'] # 132 train
Y_test_df = Y_df[Y_df.ds>'1959-12-31'] # 12 test
dataset, *_ = TimeSeriesDataset.from_df(Y_train_df)calculate_error_metrics(Y_test_df[['y']],forecasts['PatchTST'])MSE : 1530.150390625
RMSE : 39.11713678971149
MAPE : 0.07313370704650879
r2 : 0.7237732410430908
adjusted_r2 : 0.6961505651473999Y_train_df.set_index('ds',inplace =True)
forecasts.set_index('ds',inplace =True)
Y_test_df.set_index('ds',inplace =True)
plt.figure(figsize=(25, 20))<Figure size 2500x2000 with 0 Axes>y_past = Y_train_df["y"]
y_pred = forecasts['PatchTST']
y_test = Y_test_df["y"]
plt.plot(y_past, label="Past time series values")[<matplotlib.lines.Line2D object at 0x0000022F948A5CF0>]plt.plot(y_pred, label="Forecast")[<matplotlib.lines.Line2D object at 0x0000022F94882800>]plt.plot(y_test, label="Actual time series values")[<matplotlib.lines.Line2D object at 0x0000022F94883880>]plt.title('AirPassengers Forecast', fontsize=10)Text(0.5, 1.0, 'AirPassengers Forecast')plt.ylabel('Monthly Passengers', fontsize=10)Text(0, 0.5, 'Monthly Passengers')plt.xlabel('Timestamp [t]', fontsize=10)Text(0.5, 0, 'Timestamp [t]')plt.legend();
Reprogramming Large Language Model (Time-LLM)
A Large Language Model is a type of artificial intelligence model that uses machine learning techniques to process and generate human language at a scale much larger than traditional models.
These models are trained on vast amounts of text data, often encompassing entire libraries of books, website articles, social media posts, and other publicly available information. The primary function of an LLM is to predict and generate coherent and contextually relevant sequences of words, allowing it to perform tasks such as answering questions, translating languages, summarizing texts, and even generating creative writing or code.
A TIME-LLM, is a reprogramming framework to adapt large language models for time series forecasting while keeping the backbone model intact. The core idea is to reprogram the input time series into text prototype representations that are more naturally suited to language models’ capabilities.
#!pip install transformers
#!pip install psutil
#!pip install huggingface_hub
#!pip install hf_xet
#!pip install torch torchaudio torchvision --index-url https://download.pytorch.org/whl/nightly/cu128 --force-reinstall --no-cache-dir
import torch
import psutil
import platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
from neuralforecast import NeuralForecast
from neuralforecast.models import TimeLLM
from neuralforecast.utils import AirPassengersPanel, augment_calendar_df
from transformers import GPT2Config, GPT2Model, GPT2Tokenizer
use_cuda = torch.cuda.is_available()
if use_cuda:
print('__CUDNN VERSION:', torch.backends.cudnn.version())
print('__Number CUDA Devices:', torch.cuda.device_count())
print('__CUDA Device Name:',torch.cuda.get_device_name(0))
print('__CUDA Device Total Memory [GB]:',torch.cuda.get_device_properties(0).total_memory/1e9)
mem = psutil.virtual_memory()
print("Available Memory:")Available Memory:print(" Total:", mem.total / (1024 ** 2), "MB") Total: 65384.01953125 MBcpu_count = psutil.cpu_count()
cpu_count_logical = psutil.cpu_count(logical=True)
print("\nCPU Details:")
CPU Details:print(" Physical Cores:", cpu_count) Physical Cores: 32print(" Logical Cores:", cpu_count_logical) Logical Cores: 32platform_info = platform.platform()
print("\nPlatform:", platform_info)
Platform: Windows-10-10.0.26100-SP0from neuralforecast.utils import AirPassengersDF
Y_df = AirPassengersDF
Y_df = Y_df.reset_index(drop=True)
Y_df.head() unique_id ds y
0 1.0 1949-01-31 112.0
1 1.0 1949-02-28 118.0
2 1.0 1949-03-31 132.0
3 1.0 1949-04-30 129.0
4 1.0 1949-05-31 121.0Y_train_df = Y_df[Y_df.ds<='1959-12-31'] # 132 train
Y_test_df = Y_df[Y_df.ds>'1959-12-31'] # 12 test
gpt2_config = GPT2Config.from_pretrained('gpt2')
gpt2 = GPT2Model.from_pretrained('gpt2', config=gpt2_config)
gpt2_tokenizer = GPT2Tokenizer.from_pretrained('gpt2')
prompt_prefix = "The dataset contains data on monthly air passengers. There is a yearly seasonality"forecasts1.head() unique_id ds TimeLLM
0 1.0 1960-01-31 411.856049
1 1.0 1960-02-29 435.284607
2 1.0 1960-03-31 442.377167
3 1.0 1960-04-30 433.812866
4 1.0 1960-05-31 442.045441calculate_error_metrics(Y_test_df[['y']],forecasts1['TimeLLM'])MSE : 5853.29833984375
RMSE : 76.50685158757841
MAPE : 0.10379436612129211
r2 : -0.05665266513824463
adjusted_r2 : -0.1623179316520691Y_train_df.set_index('ds',inplace =True)
forecasts1.set_index('ds',inplace =True)
Y_test_df.set_index('ds',inplace =True)
plt.figure(figsize=(25, 20))<Figure size 2500x2000 with 0 Axes>y_past = Y_train_df["y"]
y_pred = forecasts1['TimeLLM']
y_test = Y_test_df["y"]
plt.plot(y_past, label="Past time series values")[<matplotlib.lines.Line2D object at 0x0000022F94984D90>]plt.plot(y_pred, label="Forecast")[<matplotlib.lines.Line2D object at 0x0000022F949BD930>]plt.plot(y_test, label="Actual time series values")[<matplotlib.lines.Line2D object at 0x0000022F949BCA00>]plt.title('AirPassengers Forecast', fontsize=10)Text(0.5, 1.0, 'AirPassengers Forecast')plt.ylabel('Monthly Passengers', fontsize=10)Text(0, 0.5, 'Monthly Passengers')plt.xlabel('Timestamp [t]', fontsize=10)Text(0.5, 0, 'Timestamp [t]')plt.legend();
A Time Series LLM for Universal Forecasting (MOIRAI)
The Masked EncOder-based UnIveRsAl Time Series Forecasting Transformer, or in short MOIRAI, is a result of novel enhancements made to the traditional time series forecasting transformer. The current version of the model was developed by training on a dataset containing observations taken across nine domains and has around 27B observations. In our experience, we found MOIRAI useful in scenarios involving zeroshot forecasting. In some use cases, it was providing results on par with multi-shot forecasting models. MOIRAI is trained and made available in three sizes, MOIRAISmall, MOIRAIBase, and MOIRAILarge, with 14m, 91m, and 311m parameters, respectively.
#!pip install gluonts
#!pip install torch
#!pip install matplotlib
#!pip install einops
#!pip install gluonts
#!pip install uni2ts
#!pip install scikit-learn