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.
- Update Oct/2017: Added a new example showing how to train on multiple prior time steps due to popular demand.
- Update Sep/2018: Updated link to dataset.
Tutorial Overview
This tutorial is divided into 3 parts; they are:
- Air Pollution Forecasting
- Basic Data Preparation
- 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:
Need help with Deep Learning for Time Series?
Take my free 7-day email crash course now (with sample code).
Click to sign-up and also get a free PDF Ebook version of the course.
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:
- No: row number
- year: year of data in this row
- month: month of data in this row
- day: day of data in this row
- hour: hour of data in this row
- pm2.5: PM2.5 concentration
- DEWP: Dew Point
- TEMP: Temperature
- PRES: Pressure
- cbwd: Combined wind direction
- Iws: Cumulated wind speed
- Is: Cumulated hours of snow
- 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.
You can download the dataset from the UCI Machine Learning Repository.
Update, I have mirrored the dataset here because UCI has become unreliable:
Download the dataset and place it in your current working directory with the filename “raw.csv“.
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):returndatetime.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()forgroup ingroups: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=1iftype(data)islist elsedata.shape[1]df=DataFrame(data)cols,names=list(),list()# input sequence (t-n, ... t-1)foriinrange(n_in,0,-1):cols.append(df.shift(i))names+=[('var%d(t-%d)'%(j+1,i))forjinrange(n_vars)]# forecast sequence (t, t+1, ... t+n)foriinrange(0,n_out):cols.append(df.shift(-i))ifi==0:names+=[('var%d(t)'%(j+1))forjinrange(n_vars)]else:names+=[('var%d(t+%d)'%(j+1,i))forjinrange(n_vars)]# put it all togetheragg=concat(cols,axis=1)agg.columns=names# drop rows with NaN valuesifdropnan:agg.dropna(inplace=True)returnagg# 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 |