Machine Learning approach - Data to Decision

Time series prediction, or forecasting, can be very similar to modelling and prediction with tabular data.
- A set of input variables, usually a single output variable.
- Ordinary machine learning methods can be applied.
The inclusion of time can be done by adding one or more delayed variables, possibly including the response.

Validation¶

As soon as time is part of a model, extra care needs to be taken in validation.
Cross-validation*, training-validation-test splits are still relevant.
However, training, validation and test sets need to follow time chronologically and cannot overlap.
Instead of traditional cross-validation, one can perform backtesting with a sliding or expanding window:

Figures from Roy Yang’s bloggpost on uber.com

Shipping, oil, interest rates and exchange rates¶

These data are public data from the Norwegian Bank, SSB, Eurostat and U.S. Energy Information Administration for the period 2000-2014 (monthly).
The data are available at ResearchGate and were part of a Master thesis by Raju Rimal.

# Read the FinalData sheet of the OilExchange.xlsx file using Pandas
import pandas as pd
# You may get a warning here, because the file contains pasted grahics
OilExchange = pd.read_excel('../../data/OilExchange.xlsx', sheet_name='FinalData') 
OilExchange.head()

OilExchange.columns

# Read the FinalCodeBook sheet of the OilExchange.xlsx file using Pandas
Explanations = pd.read_excel('../../data/OilExchange.xlsx', sheet_name='FinalCodeBook')
Explanations[['Variables','Label']]

Modelling without time¶

For starters, let us ignore time and build a simple prediction model for the exchange rate.
We will use scikit-learn’s Pipeline to combine standardisation (scaling) and linear regression and cross_val_predict to perform random K-fold cross-validation.

# Import Pipeline, StandardScaler, and LinearRegression from their respective modules in sklearn
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression

# Create a pipeline that scales the data and performs linear regression
pipe = Pipeline([('scaler', StandardScaler()), ('reg', LinearRegression())])

# Fit the pipeline with PerEURO as response and variables 3:-6 as predictors for the samples having True in the Testrain column
OilExchange_train = OilExchange.loc[OilExchange.Testrain==True,:].copy()
OilExchange_test = OilExchange.loc[OilExchange.Testrain==False,:].copy()
pipe.fit(OilExchange_train.loc[:, OilExchange_train.columns[3:-6]], \
         OilExchange_train.loc[:, 'PerEURO'])

# Predict the corresponding data for Testrain = False
PerEURO_pred = pipe.predict(OilExchange_test.loc[:, OilExchange.columns[3:-6]])

# Plot the predicted values against the actual values
import matplotlib.pyplot as plt
plt.scatter(OilExchange_test.loc[:, 'PerEURO'], PerEURO_pred)
plt.xlabel('Actual PerEURO')
plt.ylabel('Predicted PerEURO')
plt.title('Test data predictions')
plt.show()

# R2 for the test data
from sklearn.metrics import r2_score
r2_score(OilExchange_test.loc[:, 'PerEURO'], PerEURO_pred)

# Perform k-fold cross-validation with k=10
from sklearn.model_selection import cross_val_predict # NOTE: Not for time series!
PerEURO_cv = cross_val_predict(pipe, OilExchange_train.loc[:, OilExchange.columns[3:-6]], \
                               OilExchange_train.loc[:, 'PerEURO'], cv=10)

# Compute R^2 for PerEURO_cv
r2_cv = r2_score(OilExchange_train.loc[:, 'PerEURO'], PerEURO_cv)
print("Cross-validated R2: {:.3f}".format(r2_cv))

Backtesting¶

scikit-learn has a TimeSeriesSplit which creates segments for backtesting.
- Expanding window is the default.
- Sliding window can be applied by setting the right combination of parameters.
We will use scikit-learn’s cross_validate to perform the cross-validation based on the backtesting segments (cross_val_predict assumes that all observations will be test data at some point).

# Backtesting using scikit-learn
import numpy as np
from sklearn.model_selection import TimeSeriesSplit

# Some data
X = np.array([[1, 2], [3, 4], [1, 2], [3, 4], [1, 2], [3, 4]])
y = np.array([1, 2, 3, 4, 5, 6])

# Create time series cross-validation object with expanding window
tscv_expand = TimeSeriesSplit()
print(tscv_expand)
for i, (train_index, test_index) in enumerate(tscv_expand.split(X)):
    print(f"Fold {i}:")
    print(f"  Train: index={train_index}")
    print(f"  Test:  index={test_index}")

# Backtesting with sliding window
tscv_slide = TimeSeriesSplit(max_train_size=3, n_splits=3)
print(tscv_slide)
for i, (train_index, test_index) in enumerate(tscv_slide.split(X)):
    print(f"Fold {i}:")
    print(f"  Train: index={train_index}")
    print(f"  Test:  index={test_index}")

# Backtesting with expanding window in the OilExchange data
tscv_expand = TimeSeriesSplit(n_splits=10)

# The segments
max_train = []
for i, (train_index, test_index) in enumerate(tscv_expand.split(OilExchange_train.loc[:, 'PerEURO'])):
    print(f"Fold {i}:")
    print(f"  Train: index={train_index}")
    max_train.append(max(train_index))
    print(f"  Test:  index={test_index}")

# Backtesting using expanding window with data
from sklearn.model_selection import cross_validate
scores = cross_validate(pipe, OilExchange_train.loc[:, OilExchange_train.columns[3:-6]], \
                                 OilExchange_train.loc[:, 'PerEURO'], cv=tscv_expand, \
                                    scoring='r2', return_train_score=True)
scores

# Plot the backtesting results for train and test data, and under it ad the original data (PerEURO) as a subplot
plt.subplot(2,1,1)
plt.plot(scores['train_score'], label='Train')
plt.plot(scores['test_score'], label='Test')
plt.xlabel('Fold')
plt.ylabel('R$^2$')
plt.title('Backtesting results')
plt.axhline(0, color='gray', linestyle='--')
plt.ylim(-1.6,1.1)
plt.legend()
plt.subplot(2,1,2)
plt.plot(OilExchange_train.loc[:, 'PerEURO'])
for i in range(10):
    plt.axvline(x=max_train[i], color='gray', linestyle='--')
plt.xlabel('Time')
plt.ylabel('PerEURO')
plt.tight_layout()
plt.show()

Question: Does the behaviour make sense with regard to what is included in and predicted from the model?

# Backtesting with sliding window in the OilExchange data
tscv_slide = TimeSeriesSplit(max_train_size=45, n_splits=10)

# The segments
max_train = []
for i, (train_index, test_index) in enumerate(tscv_slide.split(OilExchange_train.loc[:, 'PerEURO'])):
    print(f"Fold {i}:")
    print(f"  Train: index={train_index}")
    max_train.append(max(train_index))
    print(f"  Test:  index={test_index}")

# Backtesting using sliding window with data
from sklearn.model_selection import cross_validate
scores = cross_validate(pipe, OilExchange_train.loc[:, OilExchange_train.columns[3:-4]], \
                                 OilExchange_train.loc[:, 'PerEURO'], cv=tscv_slide, \
                                    scoring='r2', return_train_score=True)
scores

# Plot the backtesting results for train and test data, and under it ad the original data (PerEURO) as a subplot
plt.subplot(2,1,1)
plt.plot(scores['train_score'], label='Train')
plt.plot(scores['test_score'], label='Test')
plt.xlabel('Fold')
plt.ylabel('R$^2$')
plt.title('Backtesting results')
plt.axhline(0, color='gray', linestyle='--')
plt.ylim(-1.6,1.1)
plt.legend()
plt.subplot(2,1,2)
plt.plot(OilExchange_train.loc[:, 'PerEURO'])
for i in range(10):
    plt.axvline(x=max_train[i], color='gray', linestyle='--')
plt.xlabel('Time')
plt.ylabel('PerEURO')
plt.tight_layout()
plt.show()

Question: Again; does the behaviour make sense with regard to what is included in and predicted from the model?

Exercise¶

Repeat the PerEuro predictions, but exchange LinearRegression with scikit-learns’s PLSRegression.
Check if the number of components in the PLS model has an effect on the explained variance ( $\text{R}^2$ ), either manually or using a GridSearchCV.

Including the response variable in the predictors¶

As long as the training and test sets are not overlapping, we can include the response as a predictor.
Adding the response lagged can be done as a single variable or several variables (i.e., several different lags).
We will later look at ARIMA-type models where time lag is the main mechanism for modelling.

# Add the Per Euro column to the OilExchange data but shifted 1 timepoint backwards (and backfill last value)
OilExchange_train['PerEURO_lag1'] = OilExchange_train.PerEURO.shift(1).bfill()
OilExchange_train.head()

# Backtesting using sliding window with data
from sklearn.model_selection import cross_validate #       Negative indexing is scary!      -->
scores = cross_validate(pipe, pd.concat([OilExchange_train.loc[:, OilExchange_train.columns[3:-7]], OilExchange_train["PerEURO_lag1"]], axis=1), \
                                 OilExchange_train.loc[:, 'PerEURO'], cv=tscv_slide, \
                                    scoring='r2', return_train_score=True)
scores

# Plot the backtesting results for train and test data, and under it ad the original data (PerEURO) as a subplot
plt.subplot(2,1,1)
plt.plot(scores['train_score'], label='Train')
plt.plot(scores['test_score'], label='Test')
plt.xlabel('Fold')
plt.ylabel('R^2')
plt.title('Backtesting results')
plt.axhline(0, color='gray', linestyle='--')
plt.ylim(-1.6,1.1)
plt.legend()
plt.subplot(2,1,2)
plt.plot(OilExchange_train.loc[:, 'PerEURO'])
for i in range(10):
    plt.axvline(x=max_train[i], color='gray', linestyle='--')
plt.xlabel('Time')
plt.ylabel('PerEURO')
plt.tight_layout()
plt.show()

Five lags¶

OilExchange_train['PerEURO_lag2'] = OilExchange_train.PerEURO.shift(2).bfill()
OilExchange_train['PerEURO_lag3'] = OilExchange_train.PerEURO.shift(3).bfill()
OilExchange_train['PerEURO_lag4'] = OilExchange_train.PerEURO.shift(4).bfill()
OilExchange_train['PerEURO_lag5'] = OilExchange_train.PerEURO.shift(5).bfill()
OilExchange_train.head()

# Backtesting using sliding window with data
from sklearn.model_selection import cross_validate #       Negative indexing is scary!      -->
scores = cross_validate(pipe, pd.concat([OilExchange_train.loc[:, OilExchange_train.columns[3:-11]], 
                                         OilExchange_train[["PerEURO_lag1","PerEURO_lag2","PerEURO_lag3","PerEURO_lag4","PerEURO_lag5"]]], axis=1),
                                         OilExchange_train.loc[:, 'PerEURO'], cv=tscv_slide,
                                         scoring='r2', return_train_score=True)
scores

# Plot the backtesting results for train and test data, and under it ad the original data (PerEURO) as a subplot
plt.subplot(2,1,1)
plt.plot(scores['train_score'], label='Train')
plt.plot(scores['test_score'], label='Test')
plt.xlabel('Fold')
plt.ylabel('R^2')
plt.title('Backtesting results')
plt.axhline(0, color='gray', linestyle='--')
plt.ylim(-1.6,1.1)
plt.legend()
plt.subplot(2,1,2)
plt.plot(OilExchange_train.loc[:, 'PerEURO'])
for i in range(10):
    plt.axvline(x=max_train[i], color='gray', linestyle='--')
plt.xlabel('Time')
plt.ylabel('PerEURO')
plt.show()

Correlation between time series¶

To get an impression of the connection between different variables, one can compute correlations, e.g., in the form of a correlation matrix.
If one expects one variable to affect another variable at a later time, correlation with a lag can be computed.
The degree of connection between two time series may also be dependent on time.
- A Sliding Window Correlation (SWC) shows local correlation in time windows.
- The window size (and possible lag) can be tuned for series of quick or slow changes.
Note: Correlation does not equal causation.
- There may not be a cause and effect, even though two phenomena show similar patterns. Beautifully illustrated by Tyler Vigen.
The concept of Autocorrelation will be covered later.

PerEURO_ExpNatGas_corr = np.corrcoef(OilExchange['PerEURO'], OilExchange['ExpNatGas'])
PerEURO_ExpNatGas_corr_lagged = np.corrcoef(OilExchange['PerEURO'][10:], OilExchange['ExpNatGas'][0:len(OilExchange['ExpNatGas'])-10])
print("Correlation between PerEURO and ExpNatGas: {:.3f}".format(PerEURO_ExpNatGas_corr[0,1]))
print("Correlation between PerEURO and ExpNatGas lagged 10 timepoints: {:.3f}".format(PerEURO_ExpNatGas_corr_lagged[0,1]))

# Use ipywidgets to create a slider for the lag
from ipywidgets import interact
def lagged_correlation(lag=0):
    x.index += lag
    corr = np.corrcoef(y[lag:], x[0:len(y)-lag])
    print("Correlation between {} and {} lagged {} timepoints: {:.3f}".format(x.name, y.name, lag, corr[0,1]))

x = OilExchange['ExpNatGas']
y = OilExchange['PerEURO']
interact(lagged_correlation, lag=(0,100,1)); # Semi-colon to suppress output

# Sliding window correlation with window size 45
PerEURO_ExpNatGas_SWC = OilExchange['PerEURO'].rolling(45, center=True).corr(OilExchange['ExpNatGas'])

# Plot PerEURO, ExpNatGas and PerEURO_ExpNatGas_SWC as subplots
def plot_SWC(center=22):
    plt.subplot(3,1,1)
    plt.plot(OilExchange['PerEURO'])
    plt.plot(range(center-22,center+22), OilExchange['PerEURO'][center-22:center+22], color="red")
    plt.ylabel('PerEURO')
    plt.xlim(0, len(OilExchange['PerEURO']))
    plt.subplot(3,1,2)
    plt.plot(OilExchange['ExpNatGas'])
    plt.plot(range(center-22,center+22), OilExchange['ExpNatGas'][center-22:center+22], color="red")
    plt.ylabel('ExpNatGas')
    plt.xlim(0, len(OilExchange['PerEURO']))
    plt.subplot(3,1,3)
    plt.plot(PerEURO_ExpNatGas_SWC)
    plt.plot(center, PerEURO_ExpNatGas_SWC[center], 'r.')
    plt.axhline(y=0, color='gray', linestyle=':')
    plt.ylim(-1,1)
    plt.xlim(0, len(OilExchange['PerEURO']))
    plt.xlabel('Time')
    plt.ylabel('SWC')
    plt.tight_layout()
    plt.show()

interact(plot_SWC, center=(22,len(OilExchange['PerEURO'])-23,1)) # Semi-colon to suppress output

Pandas’ rolling() and shifts¶

When applying Pandas’ rolling() function, the index is used for matching the data points.
Therefore, we need to shift the index of the ExpNatGas to achieve a lag.
Because of the sliding window, the two series do not need to match in length.

OE = OilExchange['ExpNatGas'].copy() # <- Remember to copy, to avoid changing the original data!
OE.index += 10
plt.plot(OilExchange['PerEURO'].rolling(45, center=True).corr(OE))
plt.xlim(0, len(OilExchange['PerEURO']))
plt.show()
OE.index

Exercise¶

Combine lag and sliding window correlation.
Use ipywidgets to control:
- window width
- lag
- selected variable to compare to PerEURO
- bonus: visualize the sliding window like above