39 KiB
39 KiB
# Import libraries
import warnings
import pandas as pd
import psutil
import sklearn
from statsmodels.tsa.seasonal import seasonal_decompose
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.neural_network import MLPRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.svm import SVR
from sklearn.ensemble import RandomForestRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
Data Preprocessing
# Load the dataset
df = pd.read_csv('data/monthly-sunspots.csv')
df.head()
| Month | Sunspots | |
|---|---|---|
| 0 | 1749-01 | 58.0 |
| 1 | 1749-02 | 62.6 |
| 2 | 1749-03 | 70.0 |
| 3 | 1749-04 | 55.7 |
| 4 | 1749-05 | 85.0 |
df.shape
(2820, 2)
# convert the date column to datetime
df['Month'] = pd.to_datetime(df['Month'], format='%Y-%m')
df.head()
| Month | Sunspots | |
|---|---|---|
| 0 | 1749-01-01 | 58.0 |
| 1 | 1749-02-01 | 62.6 |
| 2 | 1749-03-01 | 70.0 |
| 3 | 1749-04-01 | 55.7 |
| 4 | 1749-05-01 | 85.0 |
df.tail()
| Month | Sunspots | |
|---|---|---|
| 2815 | 1983-08-01 | 71.8 |
| 2816 | 1983-09-01 | 50.3 |
| 2817 | 1983-10-01 | 55.8 |
| 2818 | 1983-11-01 | 33.3 |
| 2819 | 1983-12-01 | 33.4 |
df.sort_values(by='Month', inplace=False).tail()
| Month | Sunspots | |
|---|---|---|
| 2815 | 1983-08-01 | 71.8 |
| 2816 | 1983-09-01 | 50.3 |
| 2817 | 1983-10-01 | 55.8 |
| 2818 | 1983-11-01 | 33.3 |
| 2819 | 1983-12-01 | 33.4 |
df.Sunspots.sum()
144570.0
# sum by years
df_grouped = df.groupby(pd.Grouper(key='Month', freq='Y')).sum()
df_grouped.head()
| Sunspots | |
|---|---|
| Month | |
| 1749-12-31 | 971.1 |
| 1750-12-31 | 1000.7 |
| 1751-12-31 | 571.9 |
| 1752-12-31 | 573.6 |
| 1753-12-31 | 368.3 |
# set the date column as the index
df_indexed = df.set_index('Month', inplace=False)
The data is already ordered by date. The first date is 1749-01-01 and the last is 1983-12-01. The sample is a month with its number of Sunspots registered. The data is a time series.
The sunspot number is a float, this means that is either an average or it keeps track of partial sunspots that end/start in a different month. Probably the latter since the total number of sunspots is an integer. But we can't be sure. This shouldn't matter for our model training anyway. Also, the total number is even probably because the sunspots appears in pairs. This could be useful for our model.
Sunspots are temporary phenomena on the Sun's photosphere that appear as spots darker than the surrounding areas. They are regions of reduced surface temperature caused by concentrations of magnetic field flux that inhibit convection. Sunspots usually appear in pairs of opposite magnetic polarity. Their number varies according to the approximately 11-year solar cycle.
The sunspot number is a crucial component of space weather. It is a measure of solar activity, and the number of sunspots on the solar surface changes over the course of the solar cycle. The solar cycle is a periodic change in the sun's activity and appearance. The cycle is about 11 years long on average. The solar cycle is marked by the increase and decrease of sunspots on the sun's surface. During the solar maximum, large numbers of sunspots appear, and during the solar minimum, very few sunspots appear. The solar cycle affects space weather, which can affect satellites and astronauts in space. The solar cycle is also responsible for the aurora borealis, or Northern Lights, in the Northern Hemisphere. The solar cycle was discovered in 1843 by Samuel Heinrich Schwabe, who noticed that the number of sunspots visible on the sun's surface changed over time. The solar cycle is also known as the sunspot cycle or the Schwabe cycle.
df_indexed.describe()
| Sunspots | |
|---|---|
| count | 2820.000000 |
| mean | 51.265957 |
| std | 43.448971 |
| min | 0.000000 |
| 25% | 15.700000 |
| 50% | 42.000000 |
| 75% | 74.925000 |
| max | 253.800000 |
The
df_indexed.info()
<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 2820 entries, 1749-01-01 to 1983-12-01
Data columns (total 1 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 Sunspots 2820 non-null float64
dtypes: float64(1)
memory usage: 44.1 KB
THere is no missing values in the dataset.
# print a boxplot
sns.boxplot(x=df_indexed['Sunspots'])
<Axes: xlabel='Sunspots'>
Also, no extreme values are present.
Data Visualization
Since the data is a time series, we can plot the data to see how the sunspots have changed over time. I will use the seaborn library to plot the data. Since the result is better than matplotlib.
-
Trend Component: It represents the long-term pattern or trend in the data. It shows the overall direction or tendency of the data over an extended period. The trend component indicates whether the series is increasing, decreasing, or remaining relatively stable.
-
Seasonal Component: It represents the repetitive pattern or seasonality in the data. It captures the periodic fluctuations that occur over a specific time period, such as daily, weekly, or yearly patterns.
-
Residual Component: It represents the remaining or leftover variation in the data after removing the seasonal and trend components. It consists of the random fluctuations, noise, or irregularities that cannot be explained by the seasonal or trend patterns.
# get rolling mean on a year basis
rolling_mean = df_indexed.rolling(window=12).mean().plot()
result = seasonal_decompose(df_indexed, model='additive')
sns.set(rc={'figure.figsize':(11, 4)})
ax = result.plot()
# same but in a shorted range
df_short = df_indexed['1820-01-01':'1870-01-01']
result = seasonal_decompose(df_short, model='additive')
sns.set(rc={'figure.figsize':(11, 4)})
ax = result.plot()
We can see a trend which start at 0 sunspots and increase to some number of sunspots. Then it decrease to 0 sunspots again. Graphically, this is a cycle seems to be around 10 years. We know that it is 11 years.
The seasonal component has a variance between -2.5 and 2.5.
The residual component is the noise around the trend and seasonal component.
# zooming a specific trend
df_trend = df_indexed['1832-01-01':'1844-01-01']
result = seasonal_decompose(df_trend, model='additive')
# plot with sns
sns.set(rc={'figure.figsize':(11, 4)})
ax = result.plot()
The residual component is the noise around the trend and seasonal component. It is vary between -50 and 50.
Fourier Transform
Apply the Fourier transformation and find the periodicity of the data.
The Fourier transform is a mathematical function that takes a time-based pattern as input and determines the overall cycle offset, rotation speed, and strength for every possible cycle in the given pattern. The Fourier transform is used to find the cyclic patterns hidden in the time series data.
I will calculate the periodicity on the train set to avoid test bias.
df_train = df_indexed['1749-01-01':'1919-12-01']
df_test = df_indexed['1920-01-01':'1983-12-01']
df_train.shape, df_test.shape
((2052, 1), (768, 1))
# Apply Fourier transformation
fourier_transform = np.fft.fft(df_train['Sunspots'])
fourier_transform.shape
(2052,)
fourier_transform[0:10]
array([93338.8 +0.j , 1468.1036331 -212.88095452j,
-1264.72141466-17793.04535383j, -7928.54620491 +6239.67068341j,
7122.83407343 +3535.55795089j, 555.94555365 +206.05242238j,
3791.90814354 -5119.83490291j, 1853.86932136 +4123.29141588j,
3675.5880518 +5597.02967995j, 963.21680018 -2620.91520629j])
# Compute the frequencies corresponding to the Fourier coefficients
n = len(df_train['Sunspots'])
frequencies = np.fft.fftfreq(n)
frequencies.shape
(2052,)
frequencies[0:10]
array([0. , 0.00048733, 0.00097466, 0.00146199, 0.00194932,
0.00243665, 0.00292398, 0.00341131, 0.00389864, 0.00438596])
# Find the periodicity (frequency with the highest amplitude)
periodicity = np.abs(frequencies[np.argmax(np.abs(fourier_transform))])
print('Periodicity =', periodicity)
Periodicity = 0.0
print(f"Periodicity: {1 / periodicity} time units")
Periodicity: inf time units
C:\Users\Administrator\AppData\Local\Temp\ipykernel_82340\2936873598.py:1: RuntimeWarning: divide by zero encountered in scalar divide
print(f"Periodicity: {1 / periodicity} time units")
A periodicity of "infinite" is not what we would expect.
# Plot the Fourier transform
plt.figure(figsize=(12, 4))
plt.plot(frequencies, np.abs(fourier_transform))
plt.xlabel('Frequency')
plt.ylabel('Amplitude')
plt.xlim(0, 0.02)
plt.show()
The code return a very high periodicity near 0. But this is wrong, the periodicity we want is the second peak we see. so I will recalculate the argmax but ignoring the first peak.
# Find the periodicity (frequency with the highest amplitude)
threshold = 5
max_index = np.argmax(np.abs(fourier_transform[threshold:])) + threshold
periodicity = 1 / np.abs(frequencies[max_index])
print('Periodicity =', periodicity)
Periodicity = 136.8
This means a periodicity if 136.8 months. Which is 11.4 years. This is the expected periodicity of the sunspots.
For the future engineering, I will use the periodicity of 11.4 years.
I will also calculate the periodicity of the seasonality.
# Plot the Fourier transform but this time only for the seasonality
result_season = seasonal_decompose(df_train, model='additive').seasonal
result_season.iloc[50:150].plot()
<Axes: xlabel='Month'>
n_season = len(result_season)
frequencies_season = np.fft.fftfreq(n_season)
fourier_transform_season = np.fft.fft(result_season)
plt.figure(figsize=(12, 4))
plt.plot(frequencies_season, np.abs(fourier_transform_season))
plt.xlabel('Frequency')
plt.ylabel('Amplitude')
plt.xlim(0, 0.1)
plt.show()
threshold = 0
max_index_season = np.argmax(np.abs(fourier_transform_season[threshold:])) + threshold
frequencies_season[max_index_season]
0.3333333333333333
It returns 0.333, but its clearly wrong looking at the graph, so, as a feature, I will use the seasonality of 0.083 which is 12.1 months.
Model Training
Test and Train Split
Data before 1920 as train.
# Split in the furrier transform section
df_train = df_indexed['1749-01-01':'1919-12-01']
df_test = df_indexed['1920-01-01':'1983-12-01']
df_test.head()
| Sunspots | |
|---|---|
| Month | |
| 1920-01-01 | 51.1 |
| 1920-02-01 | 53.9 |
| 1920-03-01 | 70.2 |
| 1920-04-01 | 14.8 |
| 1920-05-01 | 33.3 |
Baseline Model
As a baseline model, I will use the mean of the train data as the prediction for the test data.
Mean of the Train Data
# Create a baseline model
baseline_model = np.full(len(df_test), df_train['Sunspots'].mean())
baseline_model[0:10]
array([45.48674464, 45.48674464, 45.48674464, 45.48674464, 45.48674464,
45.48674464, 45.48674464, 45.48674464, 45.48674464, 45.48674464])
Evaluation Metrics
I will use the following metrics to evaluate the model:
- RMSE (Root Mean Squared Error): It is the square root of the average of squared differences between prediction and actual observation. It measures the standard deviation of residuals.
- MAE (Mean Absolute Error): It is the average of the absolute differences between prediction and actual observation. It measures the average magnitude of errors in a set of predictions, without considering their direction.
- R2 (R-Squared): It is the proportion of the variance in the dependent variable that is predictable from the independent variable(s). It measures how close the data are to the fitted regression line. Nearest to 1 is better.
# Calculate the RMSE (Root Mean Squared Error)
baseline_rmse = np.sqrt(mean_squared_error(df_test['Sunspots'], baseline_model))
print('Baseline RMSE:', baseline_rmse)
# Calculate the MAE (Mean Absolute Error)
baseline_mae = mean_absolute_error(df_test['Sunspots'], baseline_model)
print('Baseline MAE:', baseline_mae)
# Calculate the R2 score
baseline_r2 = r2_score(df_test['Sunspots'], baseline_model)
print('Baseline R2:', baseline_r2)
Baseline RMSE: 56.73590776250707
Baseline MAE: 43.43141345841456
Baseline R2: -0.16264654387789323
# Plot the baseline model
plt.figure(figsize=(12, 4))
plt.plot(df_train.index, df_train['Sunspots'], label='Train')
plt.plot(df_test.index, df_test['Sunspots'], label='Test')
plt.plot(df_test.index, baseline_model, label='Baseline model')
plt.legend(loc='upper left')
mae = mean_absolute_error(df_test, baseline_model)
mse = mean_squared_error(df_test, baseline_model)
r2 = r2_score(df_test, baseline_model)
plt.title(f"Model: Baseline", fontsize=16, fontweight='bold')
plt.text(0.5, 0.9, f"MAE: {mae:.2f}", transform=plt.gca().transAxes)
plt.text(0.5, 0.85, f"MSE: {mse:.2f}", transform=plt.gca().transAxes)
plt.text(0.5, 0.8, f"R^2: {r2:.2f}", transform=plt.gca().transAxes)
Text(0.5, 0.8, 'R^2: -0.16')
This clearly doesn't look good. Let's see if we can do better.
Train and Compare Models
I will create the features to train the model, the month and the year.
Features Engineering
Since the date itself is a bad feature, I need to create them, to train the model:
- The month: a number between 1 and 12 that indicates the month of the year.
- The partial: a number between 0 and 133 that indicate the month of the cycle. The cycle is 134.29 months long, so the partial is the number of months since the last cycle.
- The trend: the trend is the long-term increase or decrease in the data. It is the overall pattern of the data. It is the trend that the model will try to predict.
- The noise: the noise is the difference between the actual observation and the prediction of the previous model. It is the noise around the trend and seasonal component.
- The seasonality: the seasonality is the periodic component of the data. It is the repeating pattern within each year.
plt_seasonal_dec = seasonal_decompose(df_indexed['Sunspots'], model='additive', period=134, extrapolate_trend='freq', two_sided=True).plot()
df_train = df_indexed['1749-01-01':'1919-12-01']
df_test = df_indexed['1920-01-01':'1983-12-01']
def parse_df(df_: pd.DataFrame, trend_period: int = 137, seasonal_period: int = 12, window_size: int = 3) -> pd.DataFrame:
df_ret = pd.DataFrame(columns=[#'Month-total', # The month number since the beginning of the data (preserve total time information)
# 'Partial-trend-month', # The month number since the beginning of the cycle
# 'Partial-seasonal-month', # The month number since the beginning of the seasonality
# 'rolling-avg-3', # The rolling average, gives information about the sunspots numbers
# 'rolling_change', # The rolling difference, gives information about the direction
# 'Trend',
# 'Noise',
# 'Seasonality'
])
# df_ret['Month-total'] = df_.index.month + (df_.index.year - 1749) * 12
df_ret['Partial-trend-month'] = ((df_.index.year - 1749) * 12 + df_.index.month) % trend_period
df_ret['Partial-seasonal-month'] = ((df_.index.year - 1749) * 12 + df_.index.month) % seasonal_period
rolling_avg = df_['Sunspots'].rolling(window=3, min_periods=1).mean().to_frame()['Sunspots']
rolling_diff = df_['Sunspots'].diff(periods=1).to_frame()['Sunspots']
rolling_diff.fillna(0, inplace=True)
df_ret.index = df_.index
# to add the information about the speed of the change
# I wanted to get the rolling derivative or the gradient, but there is no function for that, however,
# I found out the percentage change function (only version 2.0.1+)
# which gives the percentage change between the current and a prior element(s). (it doesn't work)
# rolling_change = df_['Sunspots'].pcnt_change(window_size)
# I will calculate simply by doing: Dy/Dx where x is the time and y is the sunspots
# print i will use 3 month as the time step
# rolling sum of the last 3 months
rolling_der = df_['Sunspots'].diff(periods=1).to_frame()['Sunspots']
# fill the first 2 values with the first value 1
rolling_der.fillna(1, inplace=True)
# calculate the derivative
rolling_der = rolling_der/3
seasonal_decompose_result = seasonal_decompose(df_['Sunspots'], model='additive', period=134, extrapolate_trend='freq', two_sided=True)
for idx, row in df_ret.iterrows():
# df_ret.at[idx, 'rolling-avg-3'] = rolling_avg.at[idx]
df_ret.at[idx, 'rolling-diff'] = rolling_diff.at[idx]
df_ret.loc[idx, 'rolling-der'] = rolling_der.at[idx]
df_ret.loc[idx, 'Trend'] = seasonal_decompose_result.trend[idx]
df_ret.loc[idx, 'Seasonality'] = seasonal_decompose_result.seasonal[idx]
df_ret.loc[idx, 'Noise'] = seasonal_decompose_result.resid[idx]
return df_ret
X_train = parse_df(df_train)
X_test = parse_df(df_test)
# as y, use the next month sunspots
y_train = df_train['Sunspots'].shift(-1).dropna()
y_test = df_test['Sunspots'].shift(-1).dropna()
# count null
X_train.head()
| Partial-trend-month | Partial-seasonal-month | rolling-diff | rolling-der | Trend | Seasonality | Noise | |
|---|---|---|---|---|---|---|---|
| Month | |||||||
| 1749-01-01 | 1 | 1 | 0.0 | 0.333333 | 29.467287 | 25.208104 | 3.324609 |
| 1749-02-01 | 2 | 2 | 4.6 | 1.533333 | 29.579665 | 21.201296 | 11.819039 |
| 1749-03-01 | 3 | 3 | 7.4 | 2.466667 | 29.692042 | 25.258995 | 15.048963 |
| 1749-04-01 | 4 | 4 | -14.3 | -4.766667 | 29.804420 | 20.716833 | 5.178747 |
| 1749-05-01 | 5 | 5 | 29.3 | 9.766667 | 29.916798 | 20.874531 | 34.208671 |
Compare Models
def predict_next(months: int, df_X_train: pd.DataFrame(), df_y_train: pd.DataFrame() , model_: sklearn.linear_model) -> pd.DataFrame():
df_y_train = pd.DataFrame(df_y_train, columns=['Sunspots'])
# I will use the last row of the train set as the first prediction input
for i in range(months):
# this value correspond to the next month
sunspot_pred = model_.predict( df_X_train.tail(1) )
# get the last date of the df_train and add one month
last_date = df_train.tail(1).index + pd.DateOffset(months=1)
# create a df, date as index and prediction as value to the prediction set
row_pred = pd.DataFrame(data=[sunspot_pred], columns=df_y_train.columns , index=last_date)
# add the prediction row at the at to the prediction set in the only column
df_y_train = pd.concat([df_y_train, row_pred])
# also, prepend the very first sunspot that has been removed with the shift
first_sunspot = df_train['Sunspots'].head(1)
first_sunspot = pd.DataFrame(first_sunspot, columns=['Sunspots'])
df_y_train = pd.concat([first_sunspot, df_y_train])
# create the new dataframe for the prediction set based on the new prediction
df_X_train = parse_df(df_y_train)
# then shift y again
df_y_train = df_y_train.shift(-1).dropna()
return df_y_train[-months:]
import sklearn
# Models to use
models = [
LinearRegression(),
DecisionTreeRegressor(),
RandomForestRegressor(),
SVR(kernel='poly'), # SVM regressor
KNeighborsRegressor(),
MLPRegressor(),
]
months_to_predict_array = [3, 6, 12]
for months_to_predict in months_to_predict_array:
for model in models:
X_train_model = X_train.copy()
X_test_model = X_test.copy()
y_train_model = y_train.copy()
y_test_model = y_test.copy()
# fit the model, don't use the last row since has no truth
model = model.fit(X_train_model[:-1], y_train_model)
model_name = model.__class__.__name__
y_pred = predict_next(months_to_predict, X_train_model, y_train_model, model)
plt.figure(figsize=(12, 4))
# plt.plot(X_train_model[:-1].index, y_train_model, label='Train')
plt.plot(X_test_model[:months_to_predict].index, y_test_model[:months_to_predict], label='Truth')
plt.plot(X_test_model[:months_to_predict].index, y_pred, label=model_name)
plt.legend(loc='upper left')
mae = mean_absolute_error(y_test_model[:months_to_predict], y_pred)
mse = mean_squared_error(y_test_model[:months_to_predict], y_pred)
r2 = r2_score(y_test_model[:months_to_predict], y_pred)
plt.title(f"Model: {model_name}", fontsize=16, fontweight='bold')
plt.text(0.5, 0.9, f"MAE: {mae:.2f}", transform=plt.gca().transAxes)
plt.text(0.5, 0.85, f"MSE: {mse:.2f}", transform=plt.gca().transAxes)
plt.text(0.5, 0.8, f"R^2: {r2:.2f}", transform=plt.gca().transAxes)
C:\ProgramData\anaconda3\envs\MachineLearning\lib\site-packages\sklearn\neural_network\_multilayer_perceptron.py:686: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (200) reached and the optimization hasn't converged yet.
warnings.warn(
C:\ProgramData\anaconda3\envs\MachineLearning\lib\site-packages\sklearn\neural_network\_multilayer_perceptron.py:686: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (200) reached and the optimization hasn't converged yet.
warnings.warn(
C:\ProgramData\anaconda3\envs\MachineLearning\lib\site-packages\sklearn\neural_network\_multilayer_perceptron.py:686: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (200) reached and the optimization hasn't converged yet.
warnings.warn(
Conclusion
The linear regression with default value has a perfect score, so I think I could be happy with this result and stop here.
Optional : Darts
Tring using the Darts library.
from darts import TimeSeries
df_darts = pd.read_csv('data/monthly-sunspots.csv')
df_darts.head()
| Month | Sunspots | |
|---|---|---|
| 0 | 1749-01 | 58.0 |
| 1 | 1749-02 | 62.6 |
| 2 | 1749-03 | 70.0 |
| 3 | 1749-04 | 55.7 |
| 4 | 1749-05 | 85.0 |
series = TimeSeries.from_dataframe(df_darts, 'Month', 'Sunspots')
series.plot()
from darts.dataprocessing.transformers import Scaler
scaler = Scaler()
series_scaled = scaler.fit_transform(series)
series_scaled.plot()
Splitting the series
train_series, test_series = series.split_before(pd.Timestamp('19200101'))
train_series.plot(label='training')
test_series.plot(label='test')
Missing values
# example code
from darts.utils.missing_values import fill_missing_values
values = np.arange(50, step=0.5)
values[10:30] = np.nan
values[60:95] = np.nan
series_ = TimeSeries.from_values(values)
(series_ - 10).plot(label="with missing values (shifted below)")
fill_missing_values(series_).plot(label="without missing values")
Train and predict
from darts.models import NaiveSeasonal
model = NaiveSeasonal(K=1)
model.fit(train_series)
prediction = model.predict(len(test_series))
# train_series.plot(label='actual')
test_series.plot(label='test')
prediction.plot(label='forecast', lw=1)
plt.legend()
<matplotlib.legend.Legend at 0x29743e93580>
from darts.models import ExponentialSmoothing
model = ExponentialSmoothing()
model.fit(train_series)
prediction_exp = model.predict(len(test_series))
# train_series.plot(label='actual')
test_series.plot(label='actual')
prediction_exp.plot(label='forecast', lw=1)
plt.legend()
<matplotlib.legend.Legend at 0x29743a9bb80>
from darts.models import AutoARIMA
model_aarima = AutoARIMA()
model_aarima.fit(train_series)
prediction_aarima = model_aarima.predict(len(test_series))
# train_series.plot(label='actual')
test_series.plot(label='actual')
prediction_aarima.plot(label='forecast', lw=1)
plt.legend()
<matplotlib.legend.Legend at 0x29751033c10>
from darts.models import Prophet
model_prophet = Prophet()
model_prophet.fit(train_series)
prediction_prophet = model_prophet.predict(len(test_series))
# train_series.plot(label='actual')
test_series.plot(label='actual')
prediction_prophet.plot(label='forecast', lw=1)
plt.legend()
08:29:59 - cmdstanpy - INFO - Chain [1] start processing
08:29:59 - cmdstanpy - INFO - Chain [1] done processing
<matplotlib.legend.Legend at 0x29750c80670>
from darts.models import RNNModel, TCNModel, TransformerModel
model_rnn = RNNModel( input_chunk_length=12, output_chunk_length=1, n_epochs=30, model_name='LSTM')
model_rnn.fit(train_series, verbose=True)
prediction_rnn = model_rnn.predict(len(test_series))
# train_series.plot(label='actual')
test_series.plot(label='actual')
prediction_rnn.plot(label='forecast', lw=1)
plt.legend()
ignoring user defined `output_chunk_length`. RNNModel uses a fixed `output_chunk_length=1`.
GPU available: False, used: False
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
| Name | Type | Params
---------------------------------------------------
0 | criterion | MSELoss | 0
1 | train_metrics | MetricCollection | 0
2 | val_metrics | MetricCollection | 0
3 | rnn | RNN | 700
4 | V | Linear | 26
---------------------------------------------------
726 Trainable params
0 Non-trainable params
726 Total params
0.003 Total estimated model params size (MB)
Training: 0it [00:00, ?it/s]
`Trainer.fit` stopped: `max_epochs=30` reached.
GPU available: False, used: False
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
Predicting: 0it [00:00, ?it/s]
<matplotlib.legend.Legend at 0x2974e9e0430>
model_tcn = TCNModel(input_chunk_length=12, output_chunk_length=1, n_epochs=30)
model_tcn.fit(train_series, verbose=True)
prediction_tcn = model_tcn.predict(len(test_series))
# train_series.plot(label='actual')
test_series.plot(label='actual')
prediction_tcn.plot(label='forecast', lw=1)
plt.legend()
GPU available: False, used: False
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
| Name | Type | Params
----------------------------------------------------
0 | criterion | MSELoss | 0
1 | train_metrics | MetricCollection | 0
2 | val_metrics | MetricCollection | 0
3 | dropout | MonteCarloDropout | 0
4 | res_blocks | ModuleList | 92
----------------------------------------------------
92 Trainable params
0 Non-trainable params
92 Total params
0.000 Total estimated model params size (MB)
Training: 0it [00:00, ?it/s]
`Trainer.fit` stopped: `max_epochs=30` reached.
GPU available: False, used: False
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
Predicting: 0it [00:00, ?it/s]
<matplotlib.legend.Legend at 0x29751033bb0>
model_transformer = TransformerModel(input_chunk_length=12, output_chunk_length=1, n_epochs=30)
model_transformer.fit(train_series, verbose=True)
prediction_transformer = model_transformer.predict(len(test_series))
# train_series.plot(label='actual')
test_series.plot(label='actual')
prediction_transformer.plot(label='forecast', lw=1)
plt.legend()
GPU available: False, used: False
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
| Name | Type | Params
------------------------------------------------------------
0 | criterion | MSELoss | 0
1 | train_metrics | MetricCollection | 0
2 | val_metrics | MetricCollection | 0
3 | encoder | Linear | 128
4 | positional_encoding | _PositionalEncoding | 0
5 | transformer | Transformer | 548 K
6 | decoder | Linear | 65
------------------------------------------------------------
548 K Trainable params
0 Non-trainable params
548 K Total params
2.195 Total estimated model params size (MB)
Training: 0it [00:00, ?it/s]
`Trainer.fit` stopped: `max_epochs=30` reached.
GPU available: False, used: False
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
Predicting: 0it [00:00, ?it/s]
<matplotlib.legend.Legend at 0x29750ece5f0>
Evaluation
from darts.metrics import mape, mase
print('Prophet')
print('MAPE: {:.2f}%'.format(mape(prediction_prophet, test_series)))
print('MASE: {:.2f}'.format(mase(prediction_prophet, test_series, train_series)))
print('R^2: {:.2f}'.format(r2_score(test_series.values(), prediction_prophet.values())))
print('AutoARIMA')
print('MAPE: {:.2f}%'.format(mape(prediction_aarima, test_series)))
print('MASE: {:.2f}'.format(mase(prediction_aarima, test_series, train_series)))
print('R^2: {:.2f}'.format(r2_score(test_series.values(), prediction_aarima.values())))
print('Exponential Smoothing')
print('MAPE: {:.2f}%'.format(mape(prediction_exp, test_series)))
print('MASE: {:.2f}'.format(mase(prediction_exp, test_series, train_series)))
print('R^2: {:.2f}'.format(r2_score(test_series.values(), prediction_exp.values())))
print('RNN')
print('MAPE: {:.2f}%'.format(mape(prediction_rnn, test_series)))
print('MASE: {:.2f}'.format(mase(prediction_rnn, test_series, train_series)))
print('R^2: {:.2f}'.format(r2_score(test_series.values(), prediction_rnn.values())))
print('TCN')
print('MAPE: {:.2f}%'.format(mape(prediction_tcn, test_series)))
print('MASE: {:.2f}'.format(mase(prediction_tcn, test_series, train_series)))
print('R^2: {:.2f}'.format(r2_score(test_series.values(), prediction_tcn.values())))
print('Transformer')
print('MAPE: {:.2f}%'.format(mape(prediction_transformer, test_series)))
print('MASE: {:.2f}'.format(mase(prediction_transformer, test_series, train_series)))
print('R^2: {:.2f}'.format(r2_score(test_series.values(), prediction_transformer.values())))
Prophet
MAPE: 195.49%
MASE: 4.25
R^2: -0.62
AutoARIMA
MAPE: 94.51%
MASE: 3.76
R^2: -0.16
Exponential Smoothing
MAPE: 558.50%
MASE: 15.00
R^2: -14.13
RNN
MAPE: 106.90%
MASE: 3.83
R^2: -0.23
TCN
MAPE: 207.13%
MASE: 4.29
R^2: -0.67
Transformer
MAPE: 94.76%
MASE: 3.77
R^2: -0.16