Trang chủ‎ > ‎IT‎ > ‎Data Mining‎ > ‎Time Series Analysis‎ > ‎

### Multivariate Time Series Forecasting with LSTMs in Keras

Neural networks like Long Short-Term Memory (LSTM) recurrent neural networks are able to almost seamlessly model problems with multiple input variables.

This is a great benefit in time series forecasting, where classical linear methods can be difficult to adapt to multivariate or multiple input forecasting problems.

In this tutorial, you will discover how you can develop an LSTM model for multivariate time series forecasting in the Keras deep learning library.

After completing this tutorial, you will know:

• How to transform a raw dataset into something we can use for time series forecasting.
• How to prepare data and fit an LSTM for a multivariate time series forecasting problem.
• How to make a forecast and rescale the result back into the original units.

Let’s get started.

• Updated Aug/2017: Fixed a bug where yhat was compared to obs at the previous time step when calculating the final RMSE. Thanks, Songbin Xu and David Righart.

## Tutorial Overview

This tutorial is divided into 3 parts; they are:

1. Air Pollution Forecasting
2. Basic Data Preparation
3. Multivariate LSTM Forecast Model

### Python Environment

This tutorial assumes you have a Python SciPy environment installed. You can use either Python 2 or 3 with this tutorial.

You must have Keras (2.0 or higher) installed with either the TensorFlow or Theano backend.

The tutorial also assumes you have scikit-learn, Pandas, NumPy and Matplotlib installed.

If you need help with your environment, see this post:

## 1. Air Pollution Forecasting

In this tutorial, we are going to use the Air Quality dataset.

This is a dataset that reports on the weather and the level of pollution each hour for five years at the US embassy in Beijing, China.

The data includes the date-time, the pollution called PM2.5 concentration, and the weather information including dew point, temperature, pressure, wind direction, wind speed and the cumulative number of hours of snow and rain. The complete feature list in the raw data is as follows:

1. No: row number
2. year: year of data in this row
3. month: month of data in this row
4. day: day of data in this row
5. hour: hour of data in this row
6. pm2.5: PM2.5 concentration
7. DEWP: Dew Point
8. TEMP: Temperature
9. PRES: Pressure
10. cbwd: Combined wind direction
11. Iws: Cumulated wind speed
12. Is: Cumulated hours of snow
13. Ir: Cumulated hours of rain

We can use this data and frame a forecasting problem where, given the weather conditions and pollution for prior hours, we forecast the pollution at the next hour.

This dataset can be used to frame other forecasting problems.
Do you have good ideas? Let me know in the comments below.

## 2. Basic Data Preparation

The data is not ready to use. We must prepare it first.

Below are the first few rows of the raw dataset.

 123456 No,year,month,day,hour,pm2.5,DEWP,TEMP,PRES,cbwd,Iws,Is,Ir1,2010,1,1,0,NA,-21,-11,1021,NW,1.79,0,02,2010,1,1,1,NA,-21,-12,1020,NW,4.92,0,03,2010,1,1,2,NA,-21,-11,1019,NW,6.71,0,04,2010,1,1,3,NA,-21,-14,1019,NW,9.84,0,05,2010,1,1,4,NA,-20,-12,1018,NW,12.97,0,0

The first step is to consolidate the date-time information into a single date-time so that we can use it as an index in Pandas.

A quick check reveals NA values for pm2.5 for the first 24 hours. We will, therefore, need to remove the first row of data. There are also a few scattered “NA” values later in the dataset; we can mark them with 0 values for now.

The script below loads the raw dataset and parses the date-time information as the Pandas DataFrame index. The “No” column is dropped and then clearer names are specified for each column. Finally, the NA values are replaced with “0” values and the first 24 hours are removed.

The “No” column is dropped and then clearer names are specified for each column. Finally, the NA values are replaced with “0” values and the first 24 hours are removed.

 123456789101112131415161718 from pandas import read_csvfrom datetime import datetime# load datadef parse(x): return datetime.strptime(x, '%Y %m %d %H')dataset = read_csv('raw.csv',  parse_dates = [['year', 'month', 'day', 'hour']], index_col=0, date_parser=parse)dataset.drop('No', axis=1, inplace=True)# manually specify column namesdataset.columns = ['pollution', 'dew', 'temp', 'press', 'wnd_dir', 'wnd_spd', 'snow', 'rain']dataset.index.name = 'date'# mark all NA values with 0dataset['pollution'].fillna(0, inplace=True)# drop the first 24 hoursdataset = dataset[24:]# summarize first 5 rowsprint(dataset.head(5))# save to filedataset.to_csv('pollution.csv')

Running the example prints the first 5 rows of the transformed dataset and saves the dataset to “pollution.csv“.

 1234567 pollution  dew  temp   press wnd_dir  wnd_spd  snow  raindate2010-01-02 00:00:00      129.0  -16  -4.0  1020.0      SE     1.79     0     02010-01-02 01:00:00      148.0  -15  -4.0  1020.0      SE     2.68     0     02010-01-02 02:00:00      159.0  -11  -5.0  1021.0      SE     3.57     0     02010-01-02 03:00:00      181.0   -7  -5.0  1022.0      SE     5.36     1     02010-01-02 04:00:00      138.0   -7  -5.0  1022.0      SE     6.25     2     0

Now that we have the data in an easy-to-use form, we can create a quick plot of each series and see what we have.

The code below loads the new “pollution.csv” file and plots each series as a separate subplot, except wind speed dir, which is categorical.

 12345678910111213141516 from pandas import read_csvfrom matplotlib import pyplot# load datasetdataset = read_csv('pollution.csv', header=0, index_col=0)values = dataset.values# specify columns to plotgroups = [0, 1, 2, 3, 5, 6, 7]i = 1# plot each columnpyplot.figure()for group in groups: pyplot.subplot(len(groups), 1, i) pyplot.plot(values[:, group]) pyplot.title(dataset.columns[group], y=0.5, loc='right') i += 1pyplot.show()

Running the example creates a plot with 7 subplots showing the 5 years of data for each variable.

Line Plots of Air Pollution Time Series

## 3. Multivariate LSTM Forecast Model

In this section, we will fit an LSTM to the problem.

### LSTM Data Preparation

The first step is to prepare the pollution dataset for the LSTM.

This involves framing the dataset as a supervised learning problem and normalizing the input variables.

We will frame the supervised learning problem as predicting the pollution at the current hour (t) given the pollution measurement and weather conditions at the prior time step.

This formulation is straightforward and just for this demonstration. Some alternate formulations you could explore include:

• Predict the pollution for the next hour based on the weather conditions and pollution over the last 24 hours.
• Predict the pollution for the next hour as above and given the “expected” weather conditions for the next hour.

We can transform the dataset using the series_to_supervised() function developed in the blog post:

First, the “pollution.csv” dataset is loaded. The wind speed feature is label encoded (integer encoded). This could further be one-hot encoded in the future if you are interested in exploring it.

Next, all features are normalized, then the dataset is transformed into a supervised learning problem. The weather variables for the hour to be predicted (t) are then removed.

The complete code listing is provided below.

 12345678910111213141516171819202122232425262728293031323334353637383940 # convert series to supervised learningdef series_to_supervised(data, n_in=1, n_out=1, dropnan=True): n_vars = 1 if type(data) is list else data.shape[1] df = DataFrame(data) cols, names = list(), list() # input sequence (t-n, ... t-1) for i in range(n_in, 0, -1): cols.append(df.shift(i)) names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)] # forecast sequence (t, t+1, ... t+n) for i in range(0, n_out): cols.append(df.shift(-i)) if i == 0: names += [('var%d(t)' % (j+1)) for j in range(n_vars)] else: names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)] # put it all together agg = concat(cols, axis=1) agg.columns = names # drop rows with NaN values if dropnan: agg.dropna(inplace=True) return agg # load datasetdataset = read_csv('pollution.csv', header=0, index_col=0)values = dataset.values# integer encode directionencoder = LabelEncoder()values[:,4] = encoder.fit_transform(values[:,4])# ensure all data is floatvalues = values.astype('float32')# normalize featuresscaler = MinMaxScaler(feature_range=(0, 1))scaled = scaler.fit_transform(values)# frame as supervised learningreframed = series_to_supervised(scaled, 1, 1)# drop columns we don't want to predictreframed.drop(reframed.columns[[9,10,11,12,13,14,15]], axis=1, inplace=True)print(reframed.head())

Running the example prints the first 5 rows of the transformed dataset. We can see the 8 input variables (input series) and the 1 output variable (pollution level at the current hour).

 12345678910111213 var1(t-1)  var2(t-1)  var3(t-1)  var4(t-1)  var5(t-1)  var6(t-1)  \1   0.129779   0.352941   0.245902   0.527273   0.666667   0.0022902   0.148893   0.367647   0.245902   0.527273   0.666667   0.0038113   0.159960   0.426471   0.229508   0.545454   0.666667   0.0053324   0.182093   0.485294   0.229508   0.563637   0.666667   0.0083915   0.138833   0.485294   0.229508   0.563637   0.666667   0.009912    var7(t-1)  var8(t-1)   var1(t)1   0.000000        0.0  0.1488932   0.000000        0.0  0.1599603   0.000000        0.0  0.1820934   0.037037        0.0  0.1388335   0.074074        0.0  0.109658

This data preparation is simple and there is more we could explore. Some ideas you could look at include:

• One-hot encoding wind speed.
• Making all series stationary with differencing and seasonal adjustment.
• Providing more than 1 hour of input time steps.

This last point is perhaps the most important given the use of Backpropagation through time by LSTMs when learning sequence prediction problems.

### Define and Fit Model

In this section, we will fit an LSTM on the multivariate input data.

First, we must split the prepared dataset into train and test sets. To speed up the training of the model for this demonstration, we will only fit the model on the first year of data, then evaluate it on the remaining 4 years of data. If you have time, consider exploring the inverted version of this test harness.

The example below splits the dataset into train and test sets, then splits the train and test sets into input and output variables. Finally, the inputs (X) are reshaped into the 3D format expected by LSTMs, namely [samples, timesteps, features].

 123456789101112 # split into train and test setsvalues = reframed.valuesn_train_hours = 365 * 24train = values[:n_train_hours, :]test = values[n_train_hours:, :]# split into input and outputstrain_X, train_y = train[:, :-1], train[:, -1]test_X, test_y = test[:, :-1], test[:, -1]# reshape input to be 3D [samples, timesteps, features]train_X = train_X.reshape((train_X.shape[0], 1, train_X.shape[1]))test_X = test_X.reshape((test_X.shape[0], 1, test_X.shape[1]))print(train_X.shape, train_y.shape, test_X.shape, test_y.shape)

Running this example prints the shape of the train and test input and output sets with about 9K hours of data for training and about 35K hours for testing.

 1 (8760, 1, 8) (8760,) (35039, 1, 8) (35039,)

Now we can define and fit our LSTM model.

We will define the LSTM with 50 neurons in the first hidden layer and 1 neuron in the output layer for predicting pollution. The input shape will be 1 time step with 8 features.

We will use the Mean Absolute Error (MAE) loss function and the efficient Adam version of stochastic gradient descent.

The model will be fit for 50 training epochs with a batch size of 72. Remember that the internal state of the LSTM in Keras is reset at the end of each batch, so an internal state that is a function of a number of days may be helpful (try testing this).

Finally, we keep track of both the training and test loss during training by setting the validation_data argument in the fit() function. At the end of the run both the training and test loss are plotted.

 123456789101112 # design networkmodel = Sequential()model.add(LSTM(50, input_shape=(train_X.shape[1], train_X.shape[2])))model.add(Dense(1))model.compile(loss='mae', optimizer='adam')# fit networkhistory = model.fit(train_X, train_y, epochs=50, batch_size=72, validation_data=(test_X, test_y), verbose=2, shuffle=False)# plot historypyplot.plot(history.history['loss'], label='train')pyplot.plot(history.history['val_loss'], label='test')pyplot.legend()pyplot.show()

### Evaluate Model

After the model is fit, we can forecast for the entire test dataset.

We combine the forecast with the test dataset and invert the scaling. We also invert scaling on the test dataset with the expected pollution numbers.

With forecasts and actual values in their original scale, we can then calculate an error score for the model. In this case, we calculate the Root Mean Squared Error (RMSE) that gives error in the same units as the variable itself.

 123456789101112131415 # make a predictionyhat = model.predict(test_X)test_X = test_X.reshape((test_X.shape[0], test_X.shape[2]))# invert scaling for forecastinv_yhat = concatenate((yhat, test_X[:, 1:]), axis=1)inv_yhat = scaler.inverse_transform(inv_yhat)inv_yhat = inv_yhat[:,0]# invert scaling for actualtest_y = test_y.reshape((len(test_y), 1))inv_y = concatenate((test_y, test_X[:, 1:]), axis=1)inv_y = scaler.inverse_transform(inv_y)inv_y = inv_y[:,0]# calculate RMSErmse = sqrt(mean_squared_error(inv_y, inv_yhat))print('Test RMSE: %.3f' % rmse)

### Complete Example

The complete example is listed below.

NOTE: This example assumes you have prepared the data correctly, e.g. converted the downloaded “raw.csv” to the prepared “pollution.csv“. See the first part of this tutorial.

 1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495 from math import sqrtfrom numpy import concatenatefrom matplotlib import pyplotfrom pandas import read_csvfrom pandas import DataFramefrom pandas import concatfrom sklearn.preprocessing import MinMaxScalerfrom sklearn.preprocessing import LabelEncoderfrom sklearn.metrics import mean_squared_errorfrom keras.models import Sequentialfrom keras.layers import Densefrom keras.layers import LSTM # convert series to supervised learningdef series_to_supervised(data, n_in=1, n_out=1, dropnan=True): n_vars = 1 if type(data) is list else data.shape[1] df = DataFrame(data) cols, names = list(), list() # input sequence (t-n, ... t-1) for i in range(n_in, 0, -1): cols.append(df.shift(i)) names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)] # forecast sequence (t, t+1, ... t+n) for i in range(0, n_out): cols.append(df.shift(-i)) if i == 0: names += [('var%d(t)' % (j+1)) for j in range(n_vars)] else: names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)] # put it all together agg = concat(cols, axis=1) agg.columns = names # drop rows with NaN values if dropnan: agg.dropna(inplace=True) return agg # load datasetdataset = read_csv('pollution.csv', header=0, index_col=0)values = dataset.values# integer encode directionencoder = LabelEncoder()values[:,4] = encoder.fit_transform(values[:,4])# ensure all data is floatvalues = values.astype('float32')# normalize featuresscaler = MinMaxScaler(feature_range=(0, 1))scaled = scaler.fit_transform(values)# frame as supervised learningreframed = series_to_supervised(scaled, 1, 1)# drop columns we don't want to predictreframed.drop(reframed.columns[[9,10,11,12,13,14,15]], axis=1, inplace=True)print(reframed.head()) # split into train and test setsvalues = reframed.valuesn_train_hours = 365 * 24train = values[:n_train_hours, :]test = values[n_train_hours:, :]# split into input and outputstrain_X, train_y = train[:, :-1], train[:, -1]test_X, test_y = test[:, :-1], test[:, -1]# reshape input to be 3D [samples, timesteps, features]train_X = train_X.reshape((train_X.shape[0], 1, train_X.shape[1]))test_X = test_X.reshape((test_X.shape[0], 1, test_X.shape[1]))print(train_X.shape, train_y.shape, test_X.shape, test_y.shape) # design networkmodel = Sequential()model.add(LSTM(50, input_shape=(train_X.shape[1], train_X.shape[2])))model.add(Dense(1))model.compile(loss='mae', optimizer='adam')# fit networkhistory = model.fit(train_X, train_y, epochs=50, batch_size=72, validation_data=(test_X, test_y), verbose=2, shuffle=False)# plot historypyplot.plot(history.history['loss'], label='train')pyplot.plot(history.history['val_loss'], label='test')pyplot.legend()pyplot.show() # make a predictionyhat = model.predict(test_X)test_X = test_X.reshape((test_X.shape[0], test_X.shape[2]))# invert scaling for forecastinv_yhat = concatenate((yhat, test_X[:, 1:]), axis=1)inv_yhat = scaler.inverse_transform(inv_yhat)inv_yhat = inv_yhat[:,0]# invert scaling for actualtest_y = test_y.reshape((len(test_y), 1))inv_y = concatenate((test_y, test_X[:, 1:]), axis=1)inv_y = scaler.inverse_transform(inv_y)inv_y = inv_y[:,0]# calculate RMSErmse = sqrt(mean_squared_error(inv_y, inv_yhat))print('Test RMSE: %.3f' % rmse)

Running the example first creates a plot showing the train and test loss during training.

Interestingly, we can see that test loss drops below training loss. The model may be overfitting the training data. Measuring and plotting RMSE during training may shed more light on this.

Line Plot of Train and Test Loss from the Multivariate LSTM During Training

The Train and test loss are printed at the end of each training epoch. At the end of the run, the final RMSE of the model on the test dataset is printed.

We can see that the model achieves a respectable RMSE of 26.496, which is lower than an RMSE of 30 found with a persistence model.

 123456789101112 ...Epoch 46/500s - loss: 0.0143 - val_loss: 0.0133Epoch 47/500s - loss: 0.0143 - val_loss: 0.0133Epoch 48/500s - loss: 0.0144 - val_loss: 0.0133Epoch 49/500s - loss: 0.0143 - val_loss: 0.0133Epoch 50/500s - loss: 0.0144 - val_loss: 0.0133Test RMSE: 26.496

This model is not tuned. Can you do better?
Let me know your problem framing, model configuration, and RMSE in the comments below.

This section provides more resources on the topic if you are looking go deeper.

## Summary

In this tutorial, you discovered how to fit an LSTM to a multivariate time series forecasting problem.

Specifically, you learned:

• How to transform a raw dataset into something we can use for time series forecasting.
• How to prepare data and fit an LSTM for a multivariate time series forecasting problem.
• How to make a forecast and rescale the result back into the original units.

Do you have any questions?