# Understanding and implementing LSTM for multivariate time series forcasting

I want to implement an LSTM neural network to do forecasting in a multivariate time series setting and I am kind of lost when it comes to certain details regarding the structure of the network and the corresponding implementation in Python's Keras library.

I have a $d$-dimensional time series $X_t = (X^1_t,...,X^d_t)$. I want to use the last $l$ realizations to predict the next $n$ ones. To do so I use the following architecture for training (it should be $x_{t-l}$ instead of $x_{t-(l+1)}$ and $x_{t-l-1}$ instead of $x_{t-l}$):

We can see that in total we use $d+n$ time steps - $d$ in the past and $n$ in the present and future. The equations inside the LSTM cells are given by

$$f_t = \sigma(W_f[h_{t-1}, x_t]+b_f) \\ i_t = \sigma(W_i[h_{t-1}, x_t]+b_i)\\ \tilde C_t = \tanh(W_C[h_{t-1}, x_t]+b_C)\\ C_t = f_t\star C_{t-1} + i_t\star \tilde C_t\\ o_t = \sigma(W_o[h_{t-1}, x_t]+b_o)\\ h_t = o_t\star \tanh(C_t)\\$$ and the loss function for one training sample $(x_{t-l},...,x_{t+n-1})$ is given by $$\mathcal{L}(W_f, W_i, W_C, W_o, b_f, b_i, b_C, b_o) = \frac{1}{2}\sum_{k=t-(l+1)}^{t+n-1}\Vert h_{k}-X_{k+1}\Vert_{L^2(\mathbb{R}^d)}^2$$

In order for the loss function to make sense it must hold that $h_k\in\mathbb{R}^d$ and hence $W_f, W_i, W_C, W_0\in\mathbb{R}^{d\times 2d}$. Now we minimize the loss function over all training samples using Backpropagation through time.

For inference mode we then proceed as follows: We only feed in $x_{t-l},...,x_{t-1}$ for the past values and instead of the values $x_t,...,x_{t+n-1}$ to be predicted we insert $h_{t-1},...,h_{t+n-2}$:

So regarding all this I have the following questions: 1. Is this the right approach to tackle the problem? 2. How can I implement this in Keras? In keras you can enter the output dimension of the LSTM. What is this dimension? Shouldn't it be always determined by the dimension of the inputs just as explained above (s.t. the loss function makes sense). I have the feeling that keras only outputs one value at the last cell which hast the dimension entered. But then the loss function would also be different. Am I right there?