If you like Skforecast , help us giving a star on GitHub! ⭐️
More about forecasting
A time series is a succession of chronologically ordered data spaced at equal or unequal intervals. The forecasting process consists of predicting the future value of a time series, either by modeling the series solely based on its past behavior (autoregressive) or by using other external variables.
When creating a forecaster model, historical data are used to get a mathematical representation capable of predicting future values. This idea is based on a very important assumption: the future behavior of a phenomenon can be explained from its past behavior. However, this rarely happens in reality, or at least not in its entirety. For more on this, see the following definition:
$$Forecast = pattern + unexplained\;variance$$The first term of the equation refers to everything that has a repetitive character over time (trend, seasonality, cyclical factors...). The second term represents everything that influences the response variable but is not captured (explained) by the past of the time series.
The greater the importance of the first term relative to the second, the greater the probability of success when attempting to create autoregressive forecasting models. As the second term gains weight, it becomes necessary to incorporate additional variables (if any) into the model to help explain the observed behavior.
A good study of the modeled phenomenon and the ability to recognize to what extent its behavior can be explained by its past can save a lot of unnecessary effort.
This document shows an example of how to identify situations where the autoregressive forecasting process cannot achieve useful results. As an example, an attempt to predict the daily closing price of Bitcoin using machine learning methods is made. For this purpose, Skforecast is used, a simple Python library that allows, among other things, to adapt any Scikit-learn regressor to forecasting problems.
Bitcoin (₿) is a decentralized cryptocurrency that can be sent from one user to another through the bitcoin peer-to-peer network without intermediaries. Transactions are verified and recorded in a public distributed ledger called Blockchain. Bitcoins are created as a reward for a process known as mining and can be exchanged for other currencies, products, and services.
Although there may be different opinions about Bitcoin, whether as a high-risk speculative asset or, on the other hand, as a store of value, it is undeniable that it has become one of the most valuable financial assets globally. The website Infinite Market Cap shows a list of all financial assets ranked by market capitalization. Bitcoin, at the time of writing, is in the top 10. It is close to world-renowned companies such as Tesla or globally accepted safe-haven assets such as silver. The growing interest in Bitcoin, and the world of cryptocurrencies, makes it an interesting phenomenon to model.
The aim is to generate a forecasting model capable of predicting the price of Bitcoin. A time series is available with the opening (Open), closing (Close), maximum (High), and minimum (Low) prices of Bitcoin in US dollars (USD) from 2013-04-28 to 2022-01-01.
# Data manipulation
# ==============================================================================
import pandas as pd
import numpy as np
import datetime
# Plots
# ==============================================================================
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import plotly.express as px
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
from skforecast.plot import set_dark_theme
# Modelling and Forecasting
# ==============================================================================
import skforecast
import sklearn
import lightgbm
from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.model_selection import backtesting_forecaster
from lightgbm import LGBMRegressor
from sklearn.metrics import mean_absolute_error
# warnings
# ==============================================================================
import warnings
from skforecast.exceptions import LongTrainingWarning
warnings.filterwarnings("once")
warnings.simplefilter('ignore', category=LongTrainingWarning)
color = "\033[1m\033[38;5;208m"
print(f"{color}Version skforecast: {skforecast.__version__}")
print(f"{color}Version scikit-learn: {sklearn.__version__}")
print(f"{color}Version lightgbm: {lightgbm.__version__}")
print(f"{color}Version pandas: {pd.__version__}")
print(f"{color}Version numpy: {np.__version__}")
Data used in this document are daily records of the Bitcoin price in USD from 2013-04-28 to 2022-01-01. The dataset contains the following columns:
Date
: date of the record.
Open
: the opening price, the price at which an asset, in this case, Bitcoin, trades at the beginning of the day. (USD).
High
: the maximum price of the day, the highest price reached by Bitcoin on that day, (USD).
Low
: the minimum price of the day, the lowest price reached by the Bitcoin on that day, (USD).
Close
: the closing price, the price at which Bitcoin trades at the end of the day, (USD).
Volume
: the sum of actual trades made during the day, (USD).
Market Cap
: market capitalization, the total value of all shares of a company or, in the case of Bitcoin or another cryptocurrency, of all coins in circulation, (USD).
Note: the cryptocurrency market is uninterrupted. It operates 24 hours a day, 7 days a week. However, it is not strictly necessary that the close price coincides with the open price of the next day because of the fluctuations that the value of Bitcoin, or any cryptocurrency, may undergo during the last second of the day.
# Data download
# ==============================================================================
data = pd.read_csv(
"https://raw.githubusercontent.com/JoaquinAmatRodrigo/skforecast-datasets/main/data/bitcoin.csv"
)
data['date'] = pd.to_datetime(data['date'])
data = data.set_index('date')
data = data.asfreq('D')
data
When setting a frequency with the asfreq()
method, Pandas fills the gaps that may exist in the time series with the value of Null
to ensure the indicated frequency. Therefore, it should be checked if missing values have appeared after this transformation.
print(f'Number of rows with missing values: {data.isnull().any(axis=1).mean()}')
Halving is a programmed event and is part of the design and operation of some cryptocurrencies. The miners validate the transaction blocks of the network, in this case, Bitcoin, and each time they succeed, they receive an amount of that digital currency. This amount varies from time to time.
Every time 210,000 blocks are added, the reward in the Bitcoin blockchain change occurs. This event, called halving, occurs approximately every 4 years and reduces the coins miners receive as a reward by half.
In the history of Bitcoin, there have been 3 halvings. When Bitcoin mining started, miners received 50 BTC for successfully mining a block. In 2012, it reduced this reward to 25 BTC; in 2016, it dropped to 12.5 BTC, and in 2020 to 6.25 BTC, after the third halving. Each halving has affected the price, although it has not been in the short term after it.
It is intended to use the days remaining until the next halving and its mining rewards as exogenous variables to predict the price of Bitcoin. The next halving is estimated to occur approximately in 2024, although its exact date is unknown. The remaining blocks as of 2022-01-14 from the Coinmarketcap website, 121,400, and the average number of Bitcoin network blocks mined per day, 144 (average block time ≈ 10 minutes) are used to determine it.
Note: when incorporating predicted data as an exogenous variable, their error is introduced in the forecasting model since they are predictions.
# Dict with Bitcoin halvings info
# ==============================================================================
btc_halving = {
"halving": [0, 1, 2, 3, 4],
"date": ["2009-01-03", "2012-11-28", "2016-07-09", "2020-05-11", np.nan],
"reward": [50, 25, 12.5, 6.25, 3.125],
"halving_block_number": [0, 210000, 420000, 630000, 840000],
}
# Next halving calculation
# The remaining blocks according to the coinmarketcap.com website for
# the next halving as of 2022-01-14 are taken as a starting point
# ==============================================================================
remaining_blocks = 121400
blocks_per_day = 144
days = remaining_blocks / blocks_per_day
next_halving = pd.to_datetime('2022-01-14', format='%Y-%m-%d') + datetime.timedelta(days=days)
next_halving = next_halving.replace(microsecond=0, second=0, minute=0, hour=0)
next_halving = next_halving.strftime('%Y-%m-%d')
btc_halving['date'][-1] = next_halving
print(f'The next halving will occur on approximately: {next_halving}')
# Include rewards and countdown to next halving in dataset
# ==============================================================================
data["reward"] = np.nan
data["countdown_halving"] = np.nan
for i in range(len(btc_halving["halving"]) - 1):
# Start and end date of each halving
if btc_halving["date"][i] < data.index.min().strftime("%Y-%m-%d"):
start_date = data.index.min().strftime("%Y-%m-%d")
else:
start_date = btc_halving["date"][i]
end_date = btc_halving["date"][i + 1]
mask = (data.index >= start_date) & (data.index < end_date)
# Fill column 'reward' with mining rewards
data.loc[mask, "reward"] = btc_halving["reward"][i]
# Fill column 'countdown_halving' with remaining days
time_to_next_halving = pd.to_datetime(end_date) - pd.to_datetime(start_date)
data.loc[mask, "countdown_halving"] = np.arange(time_to_next_halving.days)[::-1][
: mask.sum()
]
# Check that the data have been created correctly
# ==============================================================================
print("Second halving:", btc_halving["date"][2])
display(data.loc["2016-07-08":"2016-07-09"])
print("")
print("Third halving:", btc_halving["date"][3])
display(data.loc["2020-05-10":"2020-05-11"])
print("")
print("Next halving:", btc_halving["date"][4])
data.tail(2)
When it is necessary to generate a forecasting model, plotting the time series values could be useful. This allows identifying patterns such as trends and seasonality.
A candlestick chart is a style of financial chart used to describe price movements of a security, derivative, or currency. The thick body shows the variation between the opening and closing price for a period, while the shadows show the minimum and maximum values reached during that period.
# Interactive candlestick chart with Plotly
# ==============================================================================
candlestick = go.Candlestick(
x=data.index,
open=data.open,
close=data.close,
low=data.low,
high=data.high,
)
fig = go.Figure(data=[candlestick])
fig.update_layout(
width=800,
height=350,
title=dict(text="<b>Bitcoin/USD Chart</b>", font=dict(size=20)),
yaxis_title=dict(text="Price (USD)", font=dict(size=15)),
margin=dict(l=10, r=20, t=80, b=20),
shapes=[
dict(
x0=btc_halving["date"][2],
x1=btc_halving["date"][2],
y0=0,
y1=1,
xref="x",
yref="paper",
line_width=2,
),
dict(
x0=btc_halving["date"][3],
x1=btc_halving["date"][3],
y0=0,
y1=1,
xref="x",
yref="paper",
line_width=2,
),
dict(
x0=btc_halving["date"][4],
x1=btc_halving["date"][4],
y0=0,
y1=1,
xref="x",
yref="paper",
line_width=2,
),
],
annotations=[
dict(
x=btc_halving["date"][2],
y=1,
xref="x",
yref="paper",
showarrow=False,
xanchor="left",
text="Second halving",
),
dict(
x=btc_halving["date"][3],
y=1,
xref="x",
yref="paper",
showarrow=False,
xanchor="left",
text="Third halving",
),
dict(
x=btc_halving["date"][4],
y=1,
xref="x",
yref="paper",
showarrow=False,
xanchor="left",
text="Fourth halving",
),
],
xaxis_rangeslider_visible=False,
)
fig.show()
The distribution of the Bitcoin closing price, variable close
, is shown.
# KDE distribution of Bitcoin closing price
# ==============================================================================
set_dark_theme()
fig, ax = plt.subplots(figsize=(6, 3))
p = sns.kdeplot(data=data, x='close', linewidth=3, ax=ax)
r = sns.rugplot(data=data, x='close', ax=ax, height=-.03, clip_on=False)
p.set(xlim=(0, None))
p.set_title('KDE Bitcoin closing price', fontsize=15)
p.set_xlabel('Price (USD)');
The study data shows a majority distribution of prices below 20,000 (USD). This stage corresponds to the period between 2013 and September 2020. Since 2021, the price has been in the range of 35,000 - 67,500 (USD).
Trying to model a time series with a highly asymmetric distribution and distinct orders of magnitude over time can cause complications. One strategy is model changes (deltas) rather than direct values to minimize these problems. It could be further simplified, indicating only whether the price increases or decreases from the previous day.
Note: when a distribution is asymmetric, modeling deltas instead of price may result in a more symmetric distribution.
# Locate corresponding data for each year
# ==============================================================================
years = list(data.index.year.unique())
df_plot = pd.DataFrame()
for year in years:
year_open = data.loc[data.index.year == year, "open"].iloc[0]
year_close = data.loc[data.index.year == year, "close"].iloc[-1]
year_low = data.loc[data.index.year == year, "low"].min()
year_high = data.loc[data.index.year == year, "high"].max()
df_plot[year] = pd.Series([year_open, year_close, year_low, year_high])
df_plot = df_plot.T
df_plot = df_plot.set_axis(["open", "close", "low", "high"], axis=1)
# Calculate % change between the open and close price of the year
# ==============================================================================
df_plot["year_change"] = 100 * (df_plot["close"] - df_plot["open"]) / df_plot["open"]
df_plot.head(3)
# Create a list of dicts with the % change annotations for the plot
# ==============================================================================
annotations_list = []
max_high = df_plot["high"].max()
for year in years:
df_aux = df_plot.loc[df_plot.index == year,]
loc_x = pd.to_datetime(df_aux.index[0], format="%Y")
loc_y = df_aux["high"].values[0] / max_high + 0.05
text = "{:.1f}%".format(df_aux["year_change"].values[0])
annotation = dict(
x=loc_x,
y=loc_y,
xref="x",
yref="paper",
showarrow=False,
xanchor="center",
text=text,
)
annotations_list.append(annotation)
# Interactive candlestick chart with Plotly
# ==============================================================================
candlestick = go.Candlestick(
x=pd.to_datetime(df_plot.index, format="%Y"),
open=df_plot.open,
close=df_plot.close,
low=df_plot.low,
high=df_plot.high,
)
fig = go.Figure(data=[candlestick])
fig.update_layout(
width=800,
height=350,
title=dict(text="<b>Bitcoin/USD yearly chart</b>", font=dict(size=20)),
yaxis_title=dict(text="Price (USD)", font=dict(size=13)),
margin=dict(l=0, r=20, t=55, b=20),
xaxis_rangeslider_visible=False,
annotations=annotations_list,
)
fig.show()
# Locate corresponding data for each month
# ==============================================================================
years = list(data.index.year.unique())
df_plot = pd.DataFrame()
for year in years:
for month in range(12):
start_date = pd.to_datetime(f"{year}-{month+1}-01", format="%Y-%m-%d")
end_date = start_date + pd.offsets.MonthBegin()
mask = (data.index >= start_date) & (data.index < end_date)
if not data.loc[mask, :].empty:
month_open = data.loc[mask, "open"].iloc[0]
month_close = data.loc[mask, "close"].iloc[-1]
month_low = data.loc[mask, "low"].min()
month_high = data.loc[mask, "high"].max()
serie = pd.Series([month_open, month_close, month_low, month_high])
df_aux = pd.DataFrame(serie, columns=[f"{str(month+1).zfill(2)}-{year}"])
if df_plot.empty:
df_plot = df_aux.copy()
else:
df_plot = pd.concat([df_plot, df_aux], axis=1)
df_plot = df_plot.T
df_plot = df_plot.set_axis(["open", "close", "low", "high"], axis=1)
# Boxplot chart for annual seasonality
# ==============================================================================
df_plot["month"] = pd.to_datetime(df_plot.index, format="%m-%Y").month
# fig 1 monthly boxplot
fig1 = px.box(
df_plot.sort_values("month"),
x="month",
y="close",
)
# fig 2 line with median data for each month
df_median = pd.DataFrame(df_plot.groupby("month")["close"].median()).reset_index()
fig2 = px.line(
df_median,
x="month",
y="close",
markers=True,
)
fig = go.Figure(data=fig1.data + fig2.data)
fig.update_layout(
width=650,
height=350,
title=dict(text="<b>BTC Price per month</b>", font=dict(size=20)),
yaxis_title=dict(text="Price (USD)", font=dict(size=13)),
xaxis=dict(tickmode="linear"),
xaxis_title=dict(text="month", font=dict(size=13)),
margin=dict(l=0, r=20, t=55, b=20),
)
fig.show()
# Boxplot chart for monthly seasonality
# ==============================================================================
data["day"] = pd.Series(data.index).dt.day.values
# fig 1 dayly boxplot
fig1 = px.box(data.sort_values("day"), x="day", y="close")
# fig 2 line with median data for each day
df_median = pd.DataFrame(data.groupby("day")["close"].median()).reset_index()
fig2 = px.line(df_median, x="day", y="close", markers=True)
fig = go.Figure(data=fig1.data + fig2.data)
fig.update_layout(
width=750,
height=350,
title=dict(text="<b>BTC Price per day of the month</b>", font=dict(size=20)),
yaxis_title=dict(text="Price (USD)", font=dict(size=13)),
xaxis=dict(tickmode="linear", tickangle=0, range=[0.5, 31.5]),
xaxis_title=dict(text="day", font=dict(size=13)),
margin=dict(l=0, r=20, t=55, b=20),
)
fig.show()
# Boxplot chart for weekly seasonality
# ==============================================================================
data["weekday"] = data.index.day_of_week + 1
# fig 1 weekly boxplot
fig1 = px.box(data.sort_values("weekday"), x="weekday", y="close")
# fig 2 line with median data for each weekday
df_median = pd.DataFrame(data.groupby("weekday")["close"].median()).reset_index()
fig2 = px.line(df_median, x="weekday", y="close", markers=True)
fig = go.Figure(data=fig1.data + fig2.data)
fig.update_layout(
width=600,
height=300,
title=dict(text="<b>BTC Price per day of the week</b>", font=dict(size=20)),
yaxis_title=dict(text="Price (USD)", font=dict(size=13)),
xaxis=dict(tickmode="linear"),
xaxis_title=dict(text="weekday", font=dict(size=13)),
margin=dict(l=0, r=20, t=55, b=20),
)
fig.show()
Usually, time series with an autoregressive pattern produce a repetitive shape (trend, seasonality, cyclical factors). Regarding Bitcoin, a certain annual seasonality appears at the end and beginning of the year, with larger variations in price. No seasonality appears in the monthly and weekly intervals, with very similar distributions.
# Autocorrelation plot
# ==============================================================================
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15, 4))
axes = axes.flat
plot_acf(data.close, ax=axes[0], linewidth = 1, lags=100)
axes[0].set_ylim(-1.1, 1.1)
plot_acf(data.close, ax=axes[1], linewidth = 1, lags=100)
axes[1].set_ylim(-1.1, 1.1)
plt.show();
# Partial autocorrelation plot
# ==============================================================================
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15, 4))
axes = axes.flat
plot_pacf(data.close, ax=axes[0], lags=20, method='ywm', linewidth = 1)
axes[0].set_ylim(-1.1, 1.1)
plot_pacf(data.close, ax=axes[1], lags=730, method='ywm', linewidth = 1)
axes[1].set_ylim(-1.1, 1.1)
plt.show();
The autocorrelation plots show that lag 1 is the only one correlated with lag 0. The following lags do not exceed the significance threshold.
When generating a predictive model, it is convenient to identify a base model, or baseline, on which to compare each of the proposed iterations or models. In the case of finance, according to random walk theory, market prices behave randomly and not as a function of their time series. Thus, the best estimate of the value $t_{n+1}$ is the value $t_{n}$ plus an unpredictable change.
Although a broad time series is available, it presents periods with highly differentiated prices, as determined in the graphical exploration. Only data from the last two years are used.
# Selection of train-test dates
# ==============================================================================
start_train = '2020-01-01 00:00:00'
end_train = '2021-06-30 23:59:59'
print(f"Complete time series : {data.index.min()} --- {data.index.max()} (n={len(data)})")
print(f"Train dates : {data.loc[start_train:end_train].index.min()} --- {data.loc[start_train:end_train].index.max()} (n={len(data.loc[start_train:end_train])})")
print(f"Test dates : {data.loc[end_train:].index.min()} --- {data.loc[end_train:].index.max()} (n={len(data.loc[end_train:])})")
As mentioned above, random walk suggests that the best estimate for the $t_{n+1}$ value is the $t_{n}$ value plus an unpredictable change. Therefore, the simplest viable model uses the value of $t_{n}$ as the prediction for the $t_{n+1}$ value. It is set as the baseline model.
# Random walk simulation
# ==============================================================================
df_rw = data[['close']].copy()
df_rw['pred_close'] = df_rw['close'].shift(1)
# Test error
# ==============================================================================
y_true = df_rw.loc[end_train:, 'close']
y_pred = df_rw.loc[end_train:, 'pred_close']
metric = mean_absolute_error(y_true, y_pred)
print(f'Test error: {metric}')
display(df_rw.loc[end_train:,].head(4))
# Plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 4))
df_rw.loc[end_train:, 'close'].plot(ax=ax, linewidth=2, label='Test')
df_rw.loc[end_train:, 'pred_close'].plot(ax=ax, linewidth=2, label='Prediction')
ax.set_title('Random Walk (Test data)')
ax.set_ylabel('Price (USD)')
ax.legend();
# DataFrame models' test error
# ==============================================================================
df_errors = pd.DataFrame({
'model': 'Base - Random Walk',
'lags': 1,
'test_error': metric,
'exog_variables': False
},
index=[0]
)
df_errors
The model test error is 1260.37. The following sections intend to generate a model capable of reducing this error.
An Autoregressive Forecaster (ForecasterAutoreg) with a LightGBM regressor, an implementation of the Gradient Boosting algorithm developed by Microsoft that usually achieves excellent results, is used to model the Bitcoin price.
Previous sections show an absence of autocorrelation beyond the first lag. One way to verify this is to train several models with an increasing amount of lags and check if the error doesn't reduce. This approach can be done by backtesting, using steps = 1
(predict only the next value in the series, $t_{n+1}$) it reproduces the random walk model methodology.
# Forecasters backtest with different lags
# ==============================================================================
lags = [1, 7, 30]
metrics = []
predictions_list = []
for lag in lags:
# Create forecaster
forecaster = ForecasterAutoreg(
regressor = LGBMRegressor(random_state=123, verbose=-1),
lags = lag,
transformer_y = None
)
# Backtest test data, 1 step
metric, predictions = backtesting_forecaster(
forecaster = forecaster,
y = data.loc[start_train:, 'close'],
initial_train_size = len(data.loc[start_train:end_train, 'close']),
fixed_train_size = True,
steps = 1,
refit = True,
metric = 'mean_absolute_error',
verbose = False,
show_progress = True,
)
metrics.append(metric.at[0, 'mean_absolute_error'])
predictions_list.append(predictions)
# Plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 4))
data.loc[end_train:, 'close'].plot(ax=ax, linewidth=2, label='Test')
# Plot test predictions for different lags
for predictions, lag in zip(predictions_list, lags):
predictions = predictions.rename(columns={'pred': f'Prediction, {lag} lags'})
predictions.plot(ax=ax)
ax.set_title('Close Price vs Predictions (Test Data)')
ax.set_ylabel('Price (USD)')
ax.legend();