<jupyter_start><jupyter_text>Gaussian Processes for Timeseries ForecastingIn this notebook we run some experiments to demonstrate how we can use Gaussian Processes in the context of time series forecasting. This material is part of a talk on [Gaussian Process for Time Series Analysis](https://de.pycon.org/program/pydata-knlnbb-gaussian-process-for-time-series-analysis-dr-juan-orduz/) presented at the [PyCon DE & PyData 2019 Conference](https://de.pycon.org/) in Berlin.**References:**- [Gaussian Processes for Machine Learning](http://www.gaussianprocess.org/gpml/), by <NAME> and <NAME>.- [Gaussian Processes for Timeseries Modelling](http://www.robots.ox.ac.uk/~sjrob/Pubs/philTransA_2012.pdf), by <NAME>, <NAME>, <NAME>, <NAME>, <NAME> and <NAME>.- [Bayesian Data Analysis](http://www.stat.columbia.edu/~gelman/book/), by <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, and <NAME>.- [scikit-learn docs: 1.7. Gaussian Processes](https://scikit-learn.org/stable/modules/gaussian_process.html). In particular I recomend the example: [Gaussian process regression (GPR) on Mauna Loa CO2 data](https://scikit-learn.org/stable/auto_examples/gaussian_process/plot_gpr_co2.htmlsphx-glr-auto-examples-gaussian-process-plot-gpr-co2-py).- Blog Posts:    - [Bayesian Regression as a Gaussian Process](https://juanitorduz.github.io/reg_bayesian_regression/)    - [An Introduction to Gaussian Process Regression](https://juanitorduz.github.io/gaussian_process_reg/)  Prepare Notebook<jupyter_code>import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Example 1In this first example we consider a seasonal stationary time series.  Generate Data- Time Variable Let us define a "time" variable $t\in \mathbb{N}$.<jupyter_code># Number of samples. 
n = 1000
# Generate "time" variable. 
t = np.arange(n)

data_df = pd.DataFrame({'t' : t})<jupyter_output><empty_output><jupyter_text>- Seasonal ComponentLet us define a function to generate seasonal components (Fourier modes).<jupyter_code># Generate seasonal variables. 
def seasonal(t, amplitude, period):
    """Generate a sinusoidal curve."""
    y1 = amplitude * np.sin((2*np.pi)*t/period) 
    return y1

# Add two seasonal components. 
data_df['s1'] = data_df['t'].apply(lambda t : seasonal(t, amplitude=2, period=40))

# Define target variable. 
data_df['y1'] = data_df['s1']<jupyter_output><empty_output><jupyter_text>Let us plot this seasonal variable:<jupyter_code>plt.rcParams['figure.figsize'] = [13, 8]

fig, ax = plt.subplots()
sns.lineplot(x='t', y='s1', data=data_df, label='s1', ax=ax) 
ax.set(title='Seasonal Component', xlabel='t', ylabel='')
ax.legend(loc='lower left');<jupyter_output><empty_output><jupyter_text>- Gaussian Noise Finally, we add some noise.<jupyter_code># Set noise standard deviation. 
sigma_n = 0.3

data_df['epsilon'] = np.random.normal(loc=0, scale=sigma_n, size=n)
# Add noise to target variable. 
data_df ['y1'] = data_df ['y1'] + data_df ['epsilon']<jupyter_output><empty_output><jupyter_text>Let us plot the resulting data:<jupyter_code>fig, ax = plt.subplots()
sns.lineplot(x='t', y='y1', data=data_df, label='y1', ax=ax) 
ax.set(title='Sample Data 1', xlabel='t', ylabel='')
ax.legend(loc='lower left');<jupyter_output><empty_output><jupyter_text>Define ModelLet us now define the kernel of the gaussian process model. We include the following kernel components ([recall](http://www.gaussianprocess.org/gpml/chapters/RW4.pdf) that the sum of kernels is againi a kernel):- [`WhiteKernel`](https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.WhiteKernel.html) to account for noise. - [`ExpSineSquared`](https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.ExpSineSquared.html) to model the periodic component.We add bounds to the kernel hyper-parameters which are optimized by maximizing the log-marginal-likelihood (see [documentation](https://scikit-learn.org/stable/modules/gaussian_process.html)).<jupyter_code>from sklearn.gaussian_process.kernels import WhiteKernel, ExpSineSquared, ConstantKernel

k0 = WhiteKernel(noise_level=0.3**2, noise_level_bounds=(0.1**2, 0.5**2))

k1 = ConstantKernel(constant_value=2) * ExpSineSquared(length_scale=1.0, periodicity=40, periodicity_bounds=(35, 45))

kernel_1  = k0 + k1<jupyter_output><empty_output><jupyter_text>Next we initialize the [`GaussianProcessRegressor`](https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.GaussianProcessRegressor.htmlsklearn.gaussian_process.GaussianProcessRegressor) object.<jupyter_code>from sklearn.gaussian_process import GaussianProcessRegressor

gp1 = GaussianProcessRegressor(
    kernel=kernel_1, 
    n_restarts_optimizer=10, 
    normalize_y=True,
    alpha=0.0
)<jupyter_output><empty_output><jupyter_text>Split Data  We now prepare and split the data for the model.<jupyter_code>X = data_df['t'].values.reshape(n ,1)
y = data_df['y1'].values.reshape(n ,1)

prop_train = 0.7
n_train = round(prop_train*n)

X_train = X[:n_train]
y_train = y[:n_train]

X_test = X[n_train:]
y_test = y[n_train:]<jupyter_output><empty_output><jupyter_text>Model Fit + PredictionsLet us fit the model and generate predictions.<jupyter_code>gp1.fit(X_train, y_train)
#Generate predictions.
y_pred, y_std = gp1.predict(X, return_std=True)

data_df['y_pred'] = y_pred
data_df['y_std'] = y_std
data_df['y_pred_lwr'] = data_df['y_pred'] - data_df['y_std']
data_df['y_pred_upr'] = data_df['y_pred'] + data_df['y_std']<jupyter_output><empty_output><jupyter_text>We plot the predictions.<jupyter_code>fig, ax = plt.subplots()

ax.fill_between(
    x=data_df['t'], 
    y1=data_df['y_pred_lwr'], 
    y2=data_df['y_pred_upr'], 
    color='black', 
    alpha=0.15, 
    label='credible_interval'
)

sns.lineplot(x='t', y='y1', data=data_df, label = 'y1', ax=ax)
sns.lineplot(x='t', y='y_pred', data=data_df, label='y_pred', color='black', ax=ax)

ax.axvline(n_train, color='red', linestyle='--', label='train_test_split')
ax.set(title='Prediction Sample 1', xlabel='t', ylabel='')
ax.legend(loc='upper left');<jupyter_output><empty_output><jupyter_text>Let us compute the $R^2$ of the prediction on the test set.<jupyter_code>gp1.score(X=X_test, y=y_test)<jupyter_output><empty_output><jupyter_text>Example 2In this example we add a linear trend component. - Linear Trend Component<jupyter_code># Generate trend component. 
def linear_trend(beta, x):
    """Scale vector by a scalar."""
    trend_comp = beta * x 
    return trend_comp

data_df['tr1'] = data_df['t'].apply(lambda x : linear_trend(0.01, x))

# Add trend to target variable y_1. 
data_df['y2'] = data_df['y1'] + data_df['tr1']<jupyter_output><empty_output><jupyter_text>Let us see the trend plot:<jupyter_code>fig, ax = plt.subplots()
sns.lineplot(x='t', y='y2', data=data_df, label='y_2', ax=ax)
ax.set(title='Sample Data 2', xlabel='t', ylabel='')
ax.legend(loc='upper left');<jupyter_output><empty_output><jupyter_text>Define ModelFor this second example we add a [`RBF`](https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.RBF.html) kernel to model the trend component.<jupyter_code>from sklearn.gaussian_process.kernels import RBF

k0 = WhiteKernel(noise_level=0.3**2, noise_level_bounds=(0.1**2, 0.5**2))

k1 = ConstantKernel(constant_value=2) * ExpSineSquared(length_scale=1.0, periodicity=40, periodicity_bounds=(35, 45))

k2 = ConstantKernel(constant_value=10, constant_value_bounds=(1e-2, 1e3)) * RBF(length_scale=100.0, length_scale_bounds=(1, 1e4)) 

kernel_2  = k0 + k1 + k2
# Define GaussianProcessRegressor object. 
gp2 = GaussianProcessRegressor(
    kernel=kernel_2, 
    n_restarts_optimizer=10, 
    normalize_y=True,
    alpha=0.0
)<jupyter_output><empty_output><jupyter_text>Split Data We split the data as above.<jupyter_code>y = data_df['y2'].values.reshape(n ,1)
y_train = y[:n_train]
y_test = y[n_train:]<jupyter_output><empty_output><jupyter_text>Model Fit + Predictions<jupyter_code>gp2.fit(X_train, y_train)
# Generate predictions.
y_pred, y_std = gp2.predict(X, return_std=True)

data_df['y_pred'] = y_pred
data_df['y_std'] = y_std
data_df['y_pred_lwr'] = data_df['y_pred'] - data_df['y_std']
data_df['y_pred_upr'] = data_df['y_pred'] + data_df['y_std']<jupyter_output><empty_output><jupyter_text>We plot the predictions.<jupyter_code>fig, ax = plt.subplots()

ax.fill_between(
    x=data_df['t'], 
    y1=data_df['y_pred_lwr'], 
    y2=data_df['y_pred_upr'], 
    color='black', 
    alpha=0.15, 
    label='credible_interval'
)

sns.lineplot(x='t', y='y2', data=data_df, label = 'y2', ax=ax)
sns.lineplot(x='t', y='y_pred', data=data_df, label='y_pred', color='black', ax=ax)

ax.axvline(n_train, color='red', linestyle='--', label='train_test_split')
ax.set(title='Prediction Sample 2', xlabel='t', ylabel='')
ax.legend(loc='upper left');<jupyter_output><empty_output><jupyter_text>Let us compute the $R^2$ of the prediction on the test set for this second model.<jupyter_code>gp2.score(X=X_test, y=y_test)<jupyter_output><empty_output><jupyter_text>Example 3In this third example we add a second seasonal component. - Add Other Seasonal Component<jupyter_code># Create other seasonal component.
data_df['s2'] = data_df['t'].apply(lambda t : seasonal(t, amplitude=1, period=13.3))
# Add to y_2.
data_df['y3'] = data_df['y2'] + data_df['s2'] 
fig, ax = plt.subplots()
sns.lineplot(x='t', y='y3', data=data_df, label='y3', ax=ax) 
ax.set(title='Sample Data 3', xlabel='t', ylabel='')
ax.legend(loc='upper left');<jupyter_output><empty_output><jupyter_text>Define ModelWe add another `ExpSineSquared` kernel to the one of Example 2.<jupyter_code>k0 = WhiteKernel(noise_level=0.3**2, noise_level_bounds=(0.1**2, 0.5**2))

k1 = ConstantKernel(constant_value=2) * ExpSineSquared(length_scale=1.0, periodicity=40, periodicity_bounds=(35, 45))

k2 = ConstantKernel(constant_value=10, constant_value_bounds=(1e-2, 1e3)) * RBF(length_scale=100.0, length_scale_bounds=(1, 1e4)) 

k3 = ConstantKernel(constant_value=1) * ExpSineSquared(length_scale=1.0, periodicity=12, periodicity_bounds=(10, 15))

kernel_3  = k0 + k1 + k2 + k3
# Define GaussianProcessRegressor object. 
gp3 = GaussianProcessRegressor(
    kernel=kernel_3, 
    n_restarts_optimizer=10, 
    normalize_y=True,
    alpha=0.0
)<jupyter_output><empty_output><jupyter_text>Split Data<jupyter_code>y = data_df['y3'].values.reshape(n ,1)
y_train = y[:n_train]
y_test = y[n_train:]<jupyter_output><empty_output><jupyter_text>Model Fit + Predictions<jupyter_code>gp3.fit(X_train, y_train)
y_pred, y_std = gp3.predict(X, return_std=True)

data_df['y_pred'] = y_pred
data_df['y_std'] = y_std
data_df['y_pred_lwr'] = data_df['y_pred'] - data_df['y_std']
data_df['y_pred_upr'] = data_df['y_pred'] + data_df['y_std']<jupyter_output><empty_output><jupyter_text>We plot the predictions:<jupyter_code>fig, ax = plt.subplots()

ax.fill_between(
    x=data_df['t'], 
    y1=data_df['y_pred_lwr'], 
    y2=data_df['y_pred_upr'], 
    color='black', 
    alpha=0.15, 
    label='credible_interval'
)

sns.lineplot(x='t', y='y3', data=data_df, label = 'y3', ax=ax)
sns.lineplot(x='t', y='y_pred', data=data_df, label='y_pred', color='black', ax=ax)

ax.axvline(n_train, color='red', linestyle='--', label='train_test_split')
ax.set(title='Prediction Sample Data 3', xlabel='t', ylabel='')
ax.legend(loc='upper left');
# Compute R-squared. 
gp3.score(X=X_test, y=y_test)<jupyter_output><empty_output><jupyter_text>Example 4In this last example we consider a non-linear trend component. - Non-Linear Trend Component<jupyter_code># Generate trend component. 
def non_linear_trend(x):
    """Scale and take square root."""
    trend_comp = 0.2 * np.power(x, 1/2)
    return trend_comp

# Compute non-linear trend. 
data_df ['tr2'] = data_df['t'].apply(non_linear_trend)

# Add trend to target variable. 
data_df ['y4'] = data_df ['y3'] + data_df ['tr2']<jupyter_output><empty_output><jupyter_text>Let us see the trend plot:<jupyter_code>fig, ax = plt.subplots()
sns.lineplot(x='t', y='y4', data=data_df, ax=ax)
ax.set(title='Sample Data 4', xlabel='t', ylabel='');<jupyter_output><empty_output><jupyter_text>Define ModelInstead of an `RBF` kernel to model the trendd, we use a [`RationalQuadratic`](https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.RationalQuadratic.html) which *can seen as a scale mixture (an infinite sum) of RBF kernels with different characteristic length-scales*.<jupyter_code>from sklearn.gaussian_process.kernels import RationalQuadratic

k0 = WhiteKernel(noise_level=0.3**2, noise_level_bounds=(0.1**2, 0.5**2))

k1 = ConstantKernel(constant_value=2) * ExpSineSquared(length_scale=1.0, periodicity=40, periodicity_bounds=(35, 45))

k2 = ConstantKernel(constant_value=10) * RationalQuadratic(length_scale=500, alpha= 50.0, alpha_bounds=(1, 1e3))

k3 = ConstantKernel(constant_value=1) * ExpSineSquared(length_scale=1.0, periodicity=12, periodicity_bounds=(10, 15))

kernel_4  = k0 + k1 + k2 + k3
# Define GaussianProcessRegressor object. 
gp4 = GaussianProcessRegressor(
    kernel=kernel_4, 
    n_restarts_optimizer=10, 
    normalize_y=True,
    alpha=0.0
)<jupyter_output><empty_output><jupyter_text>Split Data<jupyter_code>y = data_df['y4'].values.reshape(n ,1)
y_train = y[:n_train]
y_test = y[n_train:]<jupyter_output><empty_output><jupyter_text>Model Fit + Predictions<jupyter_code>gp4.fit(X_train, y_train)
y_pred, y_std = gp4.predict(X, return_std=True)

data_df['y_pred'] = y_pred
data_df['y_std'] = y_std
data_df['y_pred_lwr'] = data_df['y_pred'] - data_df['y_std']
data_df['y_pred_upr'] = data_df['y_pred'] + data_df['y_std']<jupyter_output><empty_output><jupyter_text>Finally, let us plot the prediction:<jupyter_code>fig, ax = plt.subplots()

ax.fill_between(
    x=data_df['t'], 
    y1=data_df['y_pred_lwr'], 
    y2=data_df['y_pred_upr'], 
    color='black', 
    alpha=0.15, 
    label='credible_interval'
)

sns.lineplot(x='t', y='y4', data=data_df, label = 'y4', ax=ax)
sns.lineplot(x='t', y='y_pred', data=data_df, label='y_pred', color='black', ax=ax)

ax.axvline(n_train, color='red', linestyle='--', label='train_test_split')
ax.set(title='Prediction Sample Data 4', xlabel='t', ylabel='')
ax.legend(loc='upper left');
# Compute R-squared. 
gp4.score(X=X_test, y=y_test)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup