Predicting parallel time series with multiple features I am trying to predict sales for 2 departmental stores which share similar demographic properties. My goal is to make a single LSTM model to predict sales from these parallel time series having multiple features.

My input features for training would be

+----------+-------+--------------+-------+
|   Date   | Store | DayOfTheWeek | Sales |
+----------+-------+--------------+-------+
| 1/1/2019 | A     |            2 |   100 |
| 1/2/2019 | A     |            3 |   200 |
| 1/3/2019 | A     |            4 |   150 |
| 1/1/2019 | B     |            2 |   300 |
| 1/2/2019 | B     |            3 |   550 |
| 1/3/2019 | B     |            4 |  1000 |
+----------+-------+--------------+-------+


and my output for training would be 

+----------+-------+--------------+-------+
|   Date   | Store | DayOfTheWeek | Sales |
+----------+-------+--------------+-------+
| 1/4/2019 | A     |            5 |   220 |
| 1/4/2019 | B     |            5 |   700 |
+----------+-------+--------------+-------+


Problem is that LSTM takes input as 3D i.e (n_sample, n_timesteps, n_features) and I can pass a single time series for a specific store (e.g. A)

If I had a univariate mutiple time series I can reshape my input data as follows and pass it to LSTM. 

+----------+---------+---------+
|   Date   | A_Sales | B_Sales |
+----------+---------+---------+
| 1/1/2019 |     100 |     300 |
| 1/2/2019 |     200 |     550 |
| 1/3/2019 |     150 |    1000 |
+----------+---------+---------+ 


But I need to identify how can I predict parallel multivariate time series? Is there any other way to define in Input LSTM layer that there are 2 time series  with 2 features each i.e (2*2).
Here are two ways to model the problem. The first one is simpler, the second one is more complex but closer to your original statement of the problem.

Store as an input feature

You can consider the store as a feature to pass to your LSTM. With two different stores, just add a binary input feature "store" where store A is 0 and store B is 1, for example.

Then, you can implement it this way.

from keras.models import Sequential
from keras.layers import LSTM, Dense

timesteps = 20
n_features = 5
model = Sequential()
# add +1 to n_features for the store identifier
model.add(LSTM(32, input_shape=(timesteps,n_features + 1), return_sequences=True))
model.add(LSTM(32))
model.add(Dense(1,activation="relu"))
model.compile(optimizer="rmsprop", loss="mse")


One sample of your data at a given timestep will be a vector of the form (feature_1, feature_2, ..., feature_n, store).

Multivariate time-series prediction

Here we input both time series and aim to predict next values of both stores. So you have a shared-LSTM processing store separately, then concatentate both produced embeddings, and compute the predicted values.

from keras.models import Model
from keras.layers import LSTM, Dense, Concatenate, Input

timesteps = 20
n_features = 5

# Store A and B time-series inputs
a_inputs = Input(shape=(timesteps, n_features))
b_inputs = Input(shape=(timesteps, n_features))

# Stacked LSTM
lstm_1 = LSTM(32, return_sequences=True)
lstm_2 = LSTM(32)
# Stacked LSTM on Store A 
a_embedding = lstm_1(a_inputs)
a_embedding = lstm_2(a_embedding)
# Stacked LSTM on Store B
b_embedding = lstm_1(b_inputs)
b_embedding = lstm_2(b_embedding)

# Concatenate embeddings and define model
outputs = Concatenate()([a_embedding, b_embedding])
outputs = Dense(64)(outputs)
outputs = Dense(2, activation="relu")(outputs)
model = Model([a_inputs, b_inputs], outputs)
model.compile(optimizer="rmsprop", loss="mse")


Inspired by section 7.1.5 Shared weight sharing, in Deep Learning with Python by F. Chollet