Introducing ProsperNN—a Python package for forecasting with neural networks

Error correction neural network

Like LSTM and GRU, the ECNN architecture (Zimmermann, Neuneier & Grothmann, 2001b) is based on RNNs. All three treat exogenous variables as model input and use the hidden state of the previous time step for calculating a new state. Although the initial state is learned during training and represents the best fit for the training data, it is equal for all inputs and therefore does not represent the related model state for a given input. For RNNs like GRU or LSTM this misspecification is forwarded without any correction and is still present in the prediction. For long input time series the effect of the initial state will vanish, but for short time series it can deteriorate forecast accuracy. Additionally, in each period new errors are included into the hidden state by model and data errors. While LSTMs and GRUs employ gates to enhance the hidden state, ECNNs use teacher forcing, i.e., calculate the prediction error in past time steps ( $t\le T$) and hand it to the respective next hidden state (Williams & Zipser, 1989). By doing this, the error which was present in the hidden state can be mitigated for the next time step, since there is correcting feedback from the data. This can improve predictions for short-term forecasting. Note, that teacher forcing in most literature hands over the target value directly to the next hidden state, e.g., Goodfellow, Courville & Bengio (2016, p.377f), whereas ECNNs subtract the prediction first, thereby passing the prediction error. Further explanations of ECNNs can be found in literature, see for example Zimmermann, Neuneier & Grothmann (2001a) or Zimmermann, Tietz & Grothmann (2012, p. 693f).

ECNNs are similar to echo state networks (Jaeger, 2001), in that they are both based on RNNs. In contrast to ECNNs, echo state networks provide additional connections from the input to the output nodes, connections between the output nodes, possibly different non-linear functions for the hidden and output nodes, and do not distinguish between past and future time steps. The training of the networks also differs: echo state networks use conventional teacher forcing and only train the weights between the hidden state and the output, since the echo state property of long-term memory only depends on the weights between the hidden states. While ECNNs rely on the error backpropagation algorithm (Rumelhart, Hinton & Williams, 1986) for training, echo state networks can be trained by solving a linear regression task, e.g., by ordinary least squares.

The idea of teacher forcing in ECNNs is similar to FORCE learning (Sussillo & Abbott, 2009), which also aims to adjust the network weights such that the error is as small as possible after training and hence external input becomes redundant. While ECNNs use the errors as additional input and learn how to use it to correct the hidden state, FORCE learning provides the hidden states with the estimation of the target without training these connections. A major difference between the two teacher forcing algorithms is that FORCE learning aims for a massive weight update after the first epoch and keeping updates small afterwards, while teacher forcing in ECNN is updated as other weights in the network. The training of the weights between the hidden states in ECNNs corresponds to innate model training for echo state networks (Laje & Buonomano, 2013).

Hence, the architecture and training of ECNNs is similar, but not the same, to combining FORCE learning on echo state networks with innate training of the weights between the hidden states.

Theory. We focus on the use of the prediction error in the ECNN since the rest of the architecture is equivalent to a basic RNN. To calculate the prediction error, we subtract the known observation from the model output: ${\hat{\mathbf{z}}}_t={\hat{\mathbf{y}}}_t-{\mathbf{y}}_t$. Afterwards, the error is multiplied with a matrix $\mathbf{D}\in {\mathbb{R}}^{k\times n}$ and added to the calculation of the next state. During training, the matrix $\mathbf{D}$ learns how to best incorporate the prediction error of the present period to improve the predictions of the following periods. The underlying assumption is that it is possible to learn where the cause of the error in the previous state stems from and that it is possible to mitigate it for the states of the following periods. So that the error ${\epsilon }_t$ in the state $s_t$ that causes the deviation from the prediction to the target, is reduced by the teacher forcing: ${\epsilon }_t\ge {\epsilon }_t+\mathbf{D}{\hat{\mathbf{z}}}_t$. As we can only calculate the error if the actual data is available, this is only possible for past values ( $t\le T$). For the state transition and prediction equations follows:

(1) $$\begin{array}{rl}{\mathbf{s}}_{t+1}& =\{\begin{array}{cc}tanh(\mathbf{A}{\mathbf{s}}_t+\mathbf{B}{\mathbf{x}}_{t+1}+\mathbf{D}({\hat{\mathbf{y}}}_t-{\mathbf{y}}_t))\hfill & t\le T\hfill \\ tanh(\mathbf{A}{\mathbf{s}}_t+\mathbf{B}{\mathbf{x}}_{t+1})\hfill & t>T\end{array}\\ {\hat{\mathbf{z}}}_t& =\mathbf{C}\tanh ({\mathbf{s}}_t)-{\mathbf{y}}_{t}t\le T\\ {\hat{\mathbf{y}}}_t& =\mathbf{C}\tanh ({\mathbf{s}}_t)t>T\end{array}$$where ${\mathbf{s}}_t\in {\mathbb{R}}^k$ is the state vector, ${\hat{\mathbf{z}}}_t$ the prediction error and $\mathbf{A}\in {\mathbb{R}}^{k\times k}$, $\mathbf{B}\in {\mathbb{R}}^{k\times m}$ and $\mathbf{C}\in {\mathbb{R}}^{n\times k}$ matrices. These relations can be seen graphically in Fig. 1. During training, we minimize ${\hat{\mathbf{z}}}_t$. Without data to correct the state (e.g., in the future), the hidden state can deteriorate over the horizon in case the model relies too much on the correction mechanism. Hence, ECNNs are most suited for short-term forecasting.

Implementation. Like the PyTorch implementation of RNNs, we implement cells that model one time step. The ECNN is stacked together by chaining multiple ECNN cells. For ECNN, the cells are designed such that they get the hidden state ${\mathbf{s}}_{t-1}$ and observation ${\mathbf{y}}_{t-1}$ of the last time step and the input ${\mathbf{x}}_t$ of the present time step. The cell then calculates the model error ${\hat{\mathbf{z}}}_{t-1}$. The teacher forcing is employed by first applying a torch.nn.Linear layer to transform ${\hat{\mathbf{z}}}_{t-1}$ and afterwards adding it to the new state at time $t$, ${\mathbf{s}}_t$. Finally, the cell returns the model error and the new model state ${\mathbf{s}}_t$. The implementation allows cases where ${\mathbf{x}}_t$ is available for $t=1,\dots ,T+\tau$ or $t=1,\dots ,T$ (see ‘Case Study’ page 18) by switching between future_U = True and future_U = False respectively. The whole ECNN (as depicted in Fig. 1) can be initialized with the parameters n_features_U (the dimension of the exogenous variables $m$), n_state_neurons (defining the size of the hidden state $k$), past_horizon (setting the look back of the model) and forecast_horizon (setting the number of steps to predict into the future):

ecnn = ECNN(n_features_U, n_state_neurons, past_horizon,
    forecast_horizon, future_U=True)

The class ECNN is a torch.nn.Module that can be used and trained in the same way as any other module in PyTorch. To train in mini-batches of size batchsize, we pass suitable data $\mathrm{x}\in {\mathbb{R}}^{(\mathtt{p}\mathtt{a}\mathtt{s}\mathtt{t}\mathtt{\_}\mathtt{h}\mathtt{o}\mathtt{r}\mathtt{i}\mathtt{z}\mathtt{o}\mathtt{n}+\tau )\times \mathtt{b}\mathtt{a}\mathtt{t}\mathtt{c}\mathtt{h}\mathtt{s}\mathtt{i}\mathtt{z}\mathtt{e}\times m}$ and Y $\in {\mathbb{R}}^{\mathtt{p}\mathtt{a}\mathtt{s}\mathtt{t}\mathtt{\_}\mathtt{h}\mathtt{o}\mathtt{r}\mathtt{i}\mathtt{z}\mathtt{o}\mathtt{n}\times \mathtt{b}\mathtt{a}\mathtt{t}\mathtt{c}\mathtt{h}\mathtt{s}\mathtt{i}\mathtt{z}\mathtt{e}\times n}$ to the model. It will return a torch.Tensor of the shape (past_horizon + forecast_horizon, batchsize, n) with the absolute errors in the first past_horizon entries of the first dimension and the prediction in the last forecast_horizon entries:

ecnn_output = ecnn(X, Y)
absolute_error = ecnn_output[:past_horizon]
forecast = ecnn_output[past_horizon:]

Historical consistent neural network

Most RNNs assume that the exogenous input variables are known in the past and in the future. This is often not fulfilled. Vector-auto-regressive models (VAR) (Stock & Watson, 2001) address this deficiency. They use multivariate input data and learn the linear dependencies between all time series to predict the future of all series. One non-linear adaption of this idea in the neural network context is the HCNN based on Zimmermann et al. (2011), which produces accurate forecasts in various practical use cases like wind forecasting (Rockefeller et al., 2023) or predicting commodity prices for procurement (Zimmermann, Grothmann & Von Mettenheim, 2013). As is common in most of these use cases, we consider the multivariate target time series $y_t$ with no input $x_t$. HCNNs can be considered as ECNNs with a special structure for multivariate forecasting. For example, it reduces parameters by fixing the teacher forcing connection $\mathbf{D}$ and is able to learn a representation of ${\mathbf{y}}_t$ in the upper part of ${\mathbf{s}}_t$.

Theory. Again, the output consists of the prediction errors along the past horizon and the forecasts of ${\mathbf{y}}_{\mathbf{t}}$ for the future. But this time, the matrix between state and prediction is fixed and just extracts the first part of the hidden state:

(2) $$\begin{array}{rl}{\hat{\mathbf{z}}}_t& =[{\mathbf{I}}_n,{\mathbf{0}}_{n,k-n}]\cdot {\mathbf{s}}_t-{\mathbf{y}}_{t}t\le T\\ {\hat{\mathbf{y}}}_t& =[{\mathbf{I}}_n,{\mathbf{0}}_{n,k-n}]\cdot {\mathbf{s}}_{t}t>T\end{array}$$where ${\mathbf{s}}_t\in {\mathbb{R}}^k$ is the state vector with $k\ge n$, ${\mathbf{I}}_n$ an $n\times n$ identity matrix, ${\mathbf{0}}_{n,k-n}$ an $n\times k-n$ zero matrix and $[\cdot ,\cdot ]$ the column-wise concatenation of two matrices so that $[{\mathbf{I}}_n,{\mathbf{0}}_{n,k-n}]\in {\mathbb{R}}^{n\times k}$.

Again, we have two cases for the state transition (see Eq. (3)), depending on whether the prediction error is used to correct the hidden state. A state ${\mathbf{s}}_{t+1}$ for $t\le T$ is calculated from the previous state ${\mathbf{s}}_t$ corrected by the prediction error ${\hat{\mathbf{z}}}_t$.

In the future ( $t>T$), the state follows only from the previous state (Eq. (3)). In both cases a state transition matrix $\mathbf{A}\in {\mathbb{R}}^{k\times k}$ is applied afterwards. Figure 2 depicts an HCNN architecture.

(3) $${\mathbf{s}}_{t+1}=\{\begin{array}{cc}\mathbf{A}tanh({\mathbf{s}}_t-{[{\mathbf{I}}_n,{\mathbf{0}}_{n,k-n}]}^{\mathrm{\top }}\cdot {\hat{\mathbf{z}}}_t)\hfill & t\le T\hfill \\ \mathbf{A}tanh({\mathbf{s}}_t)\hfill & t>T\hfill \end{array}$$

Unlike with the ECNN, the error correction here injects the true values like in common teacher forcing literature. Therefore, the intermediate state ${\tilde{\mathbf{s}}}_t$ contains ${\mathbf{y}}_t$ in the first $n$ entries before the hyperbolic tangent is applied, as seen in the following calculation.

(4) $${\tilde{\mathbf{s}}}_t={\mathbf{s}}_t-[{\mathbf{I}}_n,{\mathbf{0}}_{n,k-n}{]}^{\mathrm{\top }}\cdot {\hat{\mathbf{z}}}_t$$

(5) $$={\mathbf{s}}_t-[{\mathbf{I}}_n,{\mathbf{0}}_{n,k-n}{]}^{\mathrm{\top }}\cdot ({\mathbf{s}}_t^{[1,\dots ,n]}-{\mathbf{y}}_t)$$

(6) $$={\mathbf{s}}_t-[({\mathbf{s}}_t^{[1,\dots ,n]}-{\mathbf{y}}_t),{\mathbf{0}}_{k-n}{]}^{\mathrm{\top }}$$

(7) $$=[{\mathbf{y}}_t,{\mathbf{s}}_t^{[n+1,\dots ,k]}{]}^{\mathrm{\top }},$$where $s_t^{[1,\dots ,n]}$ is the state vector at time $t$ with only the first $n$ entries.

Implementation. HCNNs are implemented using cells that model one time step. The cell calculates the prediction and the next state according to Eqs. (2) and (3). If $\mathbf{y}$ is given to the cell, equations for $t\le T$ are used. Otherwise, equations for $t>T$ are applied.

The matrices $[{\mathbf{I}}_n,{\mathbf{0}}_{n,k-n}]$ and $[{\mathbf{I}}_n,{\mathbf{0}}_{n,k-n}{]}^{\mathrm{\top }}$ do not require gradients in the implementation, which keeps them unadjusted during training. However, they are registered as torch.nn.Parameter’s so that they are also moved to GPU if the model is moved to GPU. The whole HCNN is initialized with n_state_neurons (the number of state neurons $k$), n_features_Y (the number of variables to be predicted $n$), the past_horizon and the forecast_horizon.

hcnn = HCNN(n_state_neurons, n_features_Y, past_horizon,
    forecast_horizon)

The output of the HCNN has the shape (past_horizon + forecast_horizon, batchsize, n) and contains the in-sample absolute prediction errors in the first past_horizon entries of the first dimension and the forecasts in the remaining entries.

hcnn_output = hcnn(Y)
absolute_error = hcnn_output[:past_horizon]
forecast = hcnn_output[past_horizon:]

Causal-retro-causal neural network

CRCNNs are bidirectional HCNNs proposed in Zimmermann, Grothmann & Tietz (2012). The term of causality refers to the characteristic of HCNNs to model outputs as a function of hidden states which only contain information up to the output time. This can be seen as a non-linear interpretation of Granger causality (Granger, 1969). Like bidirectional RNNs (Schuster & Paliwal, 1997), CRCNNs are meant to model both information going forward in time (causal) and information going backward in time (retro-causal). This can be an advantage in macroeconomic forecasting as expectations of future values often influence present decisions, for instance, the expectation of high future share values leading to an increase in stock purchases and therefore a rise in current share values (Veronesi, 2000).

Theory. A basic CRCNN consists of one HCNN model going forward in time and one HCNN model going backward. Accordingly, the state $\mathbf{s}$ keeps the information that are necessary for $\mathbf{A}$ to reproduce the dynamics in the forward time direction on the one hand. On the other hand, ${\mathbf{s}}^{\mathrm{\prime }}$ and ${\mathbf{A}}^{\mathrm{\prime }}$ reproduce the dynamics in the backward time direction. The prediction value for the target at time $t$ is the sum of the prediction values, which combines the causal and the retro-causal HCNN. The two HCNNs are only connected via the shared targets ${\hat{\mathbf{z}}}_t$ (see Fig. 3). Accordingly, the state transition and prediction equations are:

(8) $$\begin{array}{rl}{\mathbf{s}}_{\mathbf{t}+\mathbf{1}}& =\{\begin{array}{cc}\mathbf{A}\tanh ({\mathbf{s}}_t-{[{\mathbf{I}}_n,{\mathbf{0}}_{n,k-n}]}^{\mathrm{\top }}\cdot {\hat{\mathbf{z}}}_t)\hfill & t\le T\hfill \\ \mathbf{A}\tanh ({\mathbf{s}}_t)\hfill & t>T\hfill \end{array}\\ {\mathbf{s}}_{\text{t−1}}^{\mathrm{\prime }}& =\{\begin{array}{cc}{\mathbf{A}}^{\mathrm{\prime }}\tanh ({\mathbf{s}}_t^{\mathrm{\prime }}-{[{\mathbf{I}}_n,{\mathbf{0}}_{n,k-n}]}^{\mathrm{\top }}\cdot {\hat{\mathbf{z}}}_t)\hfill & t\le T\hfill \\ {\mathbf{A}}^{\mathrm{\prime }}\tanh ({\mathbf{s}}_t^{\mathrm{\prime }})\hfill & t>T\hfill \end{array}\\ {\hat{\mathbf{z}}}_t& =[{\mathbf{I}}_n,{\mathbf{0}}_{n,k-n}]({\mathbf{s}}_t+{\mathbf{s}}_t^{\mathrm{\prime }})-{\mathbf{y}}_{t}t\le T\\ {\hat{\mathbf{y}}}_t& =[{\mathbf{I}}_n,{\mathbf{0}}_{n,k-n}]({\mathbf{s}}_t+{\mathbf{s}}_t^{\mathrm{\prime }})t>T.\end{array}$$

Parameter $k\ge n$ is the dimension of both of the forward $\mathbf{s}$ and backward ${\mathbf{s}}^{\mathrm{\prime }}$ state vectors.

Implementation. The use of the prediction error to correct the hidden states for both HCNNs will lead to loops (following the arrows on the highlighted path in Fig. 3 leads back to the start) in the CRCNN architecture and the error backpropagation algorithm will not converge. This is why we unfold the loop by adding copies of the HCNNs to the CRCNN model. We proceed sequentially, so that each prediction is the sum of predictions of two successive models. The prediction error is then only used for the second of every successive pair of HCNNs, which eliminates the loops while still using the prediction error on every branch but the first one. This is depicted in Fig. 4 with the smallest possible version of three branches. Teacher forcing is not applied to the bottom branch. However, for the remaining branches it is applied from one target either to the causal or the retro-causal branch, but never both. We call these HCNN models branches of the CRCNN model, with causal and retro-causal branches alternating. All causal branches share the same matrix weights as do all retro-causal branches. The total number of causal and retro-causal branches is given by n_branches.

As a retro-causal HCNN can be achieved by using HCNN cells and passing the data in reverse direction to a generic (causal) HCNN, no extra implementation of a CRCNN cell is necessary. Because of the weight sharing, we only need two HCNN cells to build a CRCNN.

A CRCNN can be initialized in the following way:

crcnn = CRCNN(n_state_neurons, n_features_Y, past_horizon,
    forecast_horizon, n_branches)

The extra parameter compared to an HCNN is defined as n_branches, which sets the total number of branches of the CRCNN.

Training is very similar to that of a single HCNN. The only difference is that the CRCNN output contains errors and forecasts for every pair of causal and retro-causal HCNNs in the model (so the first dimension contains n_branches − 1 entries).

For forecasting, we take the output of the last causal branch as it has profited most from the correction with the prediction error.

crcnn_output = crcnn(Y)
absolute_errors = crcnn_output[:, :past_horizon]
forecast = crcnn_output[−1, past_horizon:]

Fuzzy neural network

MLPs are considered black box predictors (Beck, 2018). Fuzzy neural networks whiten this black box by imposing an interpretable structure on an MLP. To make the numerical input interpretable, it classifies it as negative, zero or positive. This step is called fuzzification. The user can then propose rules for the network to follow, which are encoded in the network weights (fuzzy inference). This creates an interpretable output like negative, zero or positive that is transformed back to a numerical prediction (defuzzification). For an overview see Fig. 5. These types of fuzzy neural networks are known as ANFIS models. Walia, Singh & Sharma (2015) published a survey about the models, where they give an overview over previous research. Furthermore, implementations in PyTorch are provided by Lenhard & Maringer (2022) and https://github.com/jfpower/anfis-pytorch. In contrast to these packages, our implementation is based on Siekmann et al. (1999), which allows experts to include their prior knowledge in the form of rules into the fuzzy inference step. With these rules, experts are able to set relationships between exogenous variables and the target variable beforehand. During training the model weights the output of the rules, which can be seen as a degree of belief in each rule. Afterwards, experts can evaluate the rules by comparing in which rules they believed and in which the model does. We implemented the model for classification and regression, but focus on the later case in the following.

Theory. The first step of the fuzzy neural network is fuzzification which translates numerical input into a verbal description. This is achieved with the help of three so called membership functions, consisting of two logistic function and a Gaussian function:

(9) $$\begin{array}{rl}m{f}_0(x)& =m{f}_{neg\mathrm{\_}normlog}(x)=\frac{1}{1+\exp (-4(\lambda x-0.5))}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\lambda <0\\ m{f}_1(x)& =m{f}_{gauss}(x)=\exp \left(-\frac{1}{2}{\left(\frac{x}{\sigma }\right)}^2\right)\\ m{f}_2(x)& =m{f}_{normlog}(x)=\frac{1}{1+\exp (-4(\lambda x-0.5))}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\lambda >0.\end{array}$$

Each input variable is passed to all three functions. The function outputs correspond to the degree of membership of the fuzzy set represented. We consider an input to be a member of “negative” if $m{f}_{neg\mathrm{\_}normlog}(x)$ has the highest value, a member of “zero” if it is $m{f}_{gauss}(x)$, and a member of “positive” if it is $m{f}_{normlog}(x)$ (see Fig. 6). If the data was preprocessed with first differences, which is often done for time series, the interpretations of the three classes are decreasing, constant and increasing, respectively. During model training the function parameters, e.g., the width $\sigma$ of the Gaussian or the slope $\lambda$ of the Normlog, are adjusted, which changes the range for which the input variables are assigned to the specific classes. During training, the model learns what range of inputs corresponds to which class. For example, a volatile stock can be seen as constant in a broader range than a slow changing stock. Because these parameters are input specific, each input $x_i$ has own membership functions with varying parameters $m{f}_{neg\mathrm{\_}normlog,i}$ etc.

In the second step, fuzzy inference, we design rules between input and output variables that the network can follow, like: “IF ${\mathbf{x}}_{t,i}=zero$ AND ${\mathbf{x}}_{t,j}=negative$ THEN ${\mathbf{y}}_{t,i}=negative$”.

Conditions are modeled using the membership functions. For instance, “IF ${\mathbf{x}}_{t,i}=zero$” can be modeled as $\mathbf{b}({\mathbf{x}}_{t,i})=m{f}_{gauss,i}({\mathbf{x}}_{t,i})$, since the membership function has an output significantly higher than zero if ${\mathbf{x}}_{t,i}$ is close to zero. OR rules are built by using two IF rules, e.g., instead of “IF ${\mathbf{x}}_{t,i}=zero$ OR ${\mathbf{x}}_{t,j}=negative$ THEN ${\mathbf{y}}_{t,i}=negative$” use “IF ${\mathbf{x}}_{t,i}=zero$ THEN ${\mathbf{y}}_{t,i}=negative$” and “IF ${\mathbf{x}}_{t,j}=negative$ THEN ${\mathbf{y}}_{t,i}=negative$”. AND rules can be modeled using the product of membership functions, e.g., for the rule above: $\mathbf{b}({\mathbf{x}}_{t,i},{\mathbf{x}}_{t,j})=m{f}_{gauss,i}({\mathbf{x}}_{t,i})\cdot m{f}_{neg\mathrm{\_}normlog,j}({\mathbf{x}}_{t,j})$.

This can be extended to the $m$ inputs with each three membership functions (see Eq. (9)) and rewritten to:

(10) $${\mathbf{b}}_k(\mathbf{x})=\prod _{i=1}^m\prod _{l=0}^{2}m{f}_{l,i}({\mathbf{x}}_i{)}^{{\mathbf{P}}_{k,3i+l}}$$where ${\mathbf{b}}_k$ represents the degree to which the condition for the $k$-th rule is fulfilled and ${\mathbf{P}}_{k,3i+l}=1$ if the $i$-th input with the $l$-th membership function is part of the rule and ${\mathbf{P}}_{k,j}=0$ else. For every new rule, a new row is added to the matrix $\mathbf{P}$. In this way, any number $r$ of rules can be created. Since neural networks are not able to perform such computations, we have to rewrite the equation into a sum with activation functions:

(11) $${\mathbf{b}}_k(\mathbf{x})=\prod _{i=1}^m\prod _{l=0}^{2}m{f}_{l,i}({\mathbf{x}}_i{)}^{{\mathbf{P}}_{k,3i+l}}=\exp \left(\sum _{i=1}^m\sum _{l=0}^2{\mathbf{P}}_{k,3i+l}\ln \left(m{f}_{l,i}({\mathbf{x}}_i)\right)\right)$$

After calculating the degrees $\mathbf{b}(\mathbf{x})$ to which the rules are fulfilled, we multiply with a matrix $\mathbf{C}\in {\mathbb{R}}^{3\times r}$. It models the consequences of the rules by connecting the entries of $\mathbf{b}(\mathbf{x})$ to three classes negative, zero and positive:

(12) $${\mathbf{x}}_{interpretable}=\mathbf{C}\cdot softmax(\mathbf{b})$$

For each rule $i$, the column $i$ of $\mathbf{C}$ models the respective consequence by setting the connection to either one or zero. For instance, if the consequence of a rule is negative, only the entry in $\mathbf{C}$ that connects ${\mathbf{b}}_i$ and the class that represents negative is one and the entries to the classes zero and positive are zero. During training the entries of $\mathbf{C}$ are adjusted and can be seen as a measure of how much the model trusts each rule given the data. The entries of ${\mathbf{x}}_{interpretable}$ in Eq. (12) show to what degree the forecast is negative, zero and positive.

Finally, the defuzzification translates the interpretable output ${\mathbf{x}}_{interpretable}$ back to a numerical value with $\mathbf{W}\in {\mathbb{R}}^{1\times 3}$:

(13) $${\mathbf{y}}_{t+1}=\mathbf{W}{\mathbf{x}}_{interpretable}$$

Implementation. For fuzzification, we have written a torch.nn.Module for each of the three membership functions defined in Eq. (9) and another one that applies them to any input variable. So for each input variable, the parameters of the three functions are trained independently. The output of the fuzzification has one dimension more than the input since each input variable results in three return values.

To initialize a Fuzzy-NN, we first have to define the membership functions and their meaning. With the proposed structure of MemberLayer users can add and replace membership functions to enable new ways of model and data interpretation.

membership_fcts = {
 "negative": NormlogMember(negative=True),
 "zero": GaussianMember(),
 "positive": NormlogMember()}
fuzzification = Fuzzification(
 n_features_input=n_input,
 membership_fcts=membership_fcts)

In fuzzy inference the RuleManager class loads the user-defined rules from a json file. They consist of conditions and consequences, the name and position of the input names and the interpretation of the membership functions. The RuleManager automatically creates the corresponding matrices $\mathbf{P}$ and $\mathbf{C}$. The matrix $\mathbf{P}$ is implemented with a convolutional layer $torch.nn.Conv1d$, since the output of the fuzzification step is two dimensional (number of input variables times the number of rules). With these rules we can initialize the torch.nn.Linear layer that performs Eq. (11). It is possible to create an arbitrary number of rules and output classes. So the interpretation of the output is not limited to low, constant and high, but allows adaption to various use cases.

rule_manager = RuleManager(
 path="rules.json",
 rule_matrix_shape=(3m, n_input, 3),
 classification_matrix_shape=(r, 3))
fuzzy_inference = FuzzyInference(
 n_features_input=m,
 n_rules=r,
 n_output_classes=3,
 n_membership_fcts=3,
 rule_matrix=rule_manager.rule_matrix,
 classification_matrix=rule_manager.classification_matrix)

In the defuzzification, first the matrix $\mathbf{C}$ is applied so that we get the interpretable output. Afterwards, a softmax function and the matrix $\mathbf{W}$ are applied.

defuzzification = Defuzzification(n_output_classes=3)

Finally, after all layers are initialized, we can create the Fuzzy Neural Network that can be trained like any other PyTorch module:

fuzzy_nn = torch.nn.Sequential(fuzzification, fuzzy_inference,
    defuzzification)

Visualizations of forecasts and uncertainties

Ensembles are frequently used to quantify uncertainty of neural networks (Lakshminarayanan, Pritzel & Blundell, 2017; Van Schaeybroeck & Vannitsem, 2016). Köpp, Von Mettenheim & Breitner (2014) present an approach using heatmaps to visualize tiny accumulations and paths taken by many models. An implementation is only available in Java (Köpp, Von Mettenheim & Breitner, 2014). We implement this visualization in Python and demonstrate it using a single output variable as an example.

Theory. To display the heatmap continuously, we interpolate the values between the discrete time steps in the forecast. Based on these continuous forecasts, we then determine the heat of every pixel in the heatmap.

Let $d\in \mathbb{N}$ be the number of forecasts, $p\in \mathbb{N}$ the number of interpolated values between two time steps and $y_{res}\in \mathbb{N}$ the resolution in the y-axis. The number of pixels in the heatmap is therefore $y_{res}\times p\tau$. Increasing either $p$ or $y_{res}$ smoothens the heatmap, but with computational costs.

We interpolate all forecasts linearly between all time steps:

(14) $${\hat{y}}_{t+j/p}^{[i]}=\{\begin{array}{cc}{y}_T+\frac{j}{p}\cdot ({\hat{y}}_{T+1}^{[i]}-y_T)\hfill & \mathrm{i}\mathrm{f}t=T\hfill \\ {\hat{y}}_t^{[i]}+\frac{j}{p}\cdot ({\hat{y}}_{t+1}^{[i]}-{\hat{y}}_t^{[i]})\hfill & \mathrm{i}\mathrm{f}t=T+1,\dots ,T+\tau \hfill \end{array}\mathrm{f}\mathrm{o}\mathrm{r}\begin{array}{c}i=1,\dots ,d\hfill \\ j=1,\dots ,p-1\hfill \end{array}$$where $j/p$ defines the intermediate time steps between periods and ${\hat{y}}_{t+j/p}^{[i]}$ represents the interpolated value for the $i$-th forecast and time step $t+j/p$.

Next, we calculate a matrix $H\in {\mathbb{R}}^{y_{res}\times (p\tau )}$ where each entry represents the heat of one pixel in the figure. There is a column for each interpolated time step $t=T,T+\frac{1}{p},\dots ,T+\frac{p\tau }{p}$. The heat is based on the difference between the value the pixel represents and the forecasted values. Therefore, we assign each row $H_u$ in the figure a value $f(u)$ on the same scale as $y$ and the forecasts. We define the bottom pixel row ( $u=1$) as the minimum of the forecasted values and the last known value $y_T$: ${\hat{y}}_{min}=min(\{{\hat{y}}_t^{[i]}|t=T,\dots ,T+\tau ;i=1,\dots ,d\}\cup \{y_T\})$ and the top row ( $u=y_{res}$) as the maximum: ${\hat{y}}_{max}=max(\{{\hat{y}}_t^{[i]}|t=T,\dots ,T+\tau ;i=1,\dots ,d\}\cup \{y_T\})$. The values of the rows in between are interpolated linearly: $f(u)={\hat{y}}_{min}+\frac{u}{y_{res}}({\hat{y}}_{max}-{\hat{y}}_{min})$.

Each forecast contributes to the total heat of the $t$-th pixel in row $u$ with the difference between the value the pixel represents $f(u)$ and the forecast ${\hat{y}}_t^{[i]}$. The difference is weighted by the Gaussian kernel $k_{Gaussian}(\delta )={\exp }^{-{(\frac{\delta }{\sigma })}^2}$ with $\sigma \in {\mathbb{R}}^{+}$. For each pixel we sum up the heat of all forecasts at the time point:

(15) $$h_{u,v}=\sum _{i=1}^d{\exp }^{-{\left(\frac{f(u)-{\hat{y}}_v^{[i]}}{\sigma }\right)}^2}\mathrm{f}\mathrm{o}\mathrm{r}\begin{array}{rl}u& =1,\dots ,y_{res}\\ v& =T,T+\frac{1}{p},\dots ,T+\frac{p\tau }{p}\end{array}$$

Since at earlier forecast steps the heat is higher due to denser forecasts, accumulations for less dense regions are hard to identify (Köpp, Von Mettenheim & Breitner, 2014). That is why we normalize each column to have a maximum heat of 1 in each step:

(16) $$h_{u,v}^{norm}=\frac{h_{u,v}}{{max}_{u=1,\dots ,y_{res}}{h}_{u,v}}\mathrm{f}\mathrm{o}\mathrm{r}v=T,T+\frac{1}{p},\dots ,T+\frac{p\tau }{p}$$

Implementation. In order to simplify creating an ensemble, we implemented a torch.nn.Module that takes a base_model and duplicates it, using different start initializations for the weights. With the parameter initializer in Ensemble the user can choose the method for weight initialization. Initializers from torch.nn.init are possible with the default value torch.nn.init.kaiming_uniform_. With n_models we can set the size of the ensemble:

ensemble = Ensemble(model=base_model, n_models=n_models,
    initializer=kaiming_uniform_)

The variable forecasts has to be a torch.Tensor with $d$ forecasts in the rows and a prediction step $\tau$ in each column. The heatmap_forecasts function uses scipy’s interp1d function for interpolation and seaborn to visualize the heatmap (Virtanen et al., 2020; Waskom, 2021). The forecasts in the heatmap are depicted sharper with sigma close to zero and smoother with bigger values. In the function call, we have to set the parameters $d,p$ and $y_{res}$ (see ‘Theory’). The number of forecasts is taken from the number of rows of forecasts.

heatmap_forecasts(forecasts, sigma, n_interpolation=p,
    y_resolution=y_res)

Sensitivity analysis

Due to the various application fields, there is growing interest in the interpretability of neural networks for time series analysis (Rojat et al., 2021). Methods can be divided into post-hoc and ante-hoc as well as local and global ones (Rojat et al., 2021): While post-hoc methods examine the relationship of model input and output after predictions, ante-hoc methods are embedded within the model architecture, so that interpretability is included in model training. Examples for the latter are N-BEATs (Oreshkin et al., 2019) which can split the forecast into components, attention layers as applied within Transformer models (Vaswani et al., 2017), and fuzzy neural networks as discussed in ‘Fuzzy neural network’. Local interpretability approaches act on single samples, e.g., a single forecast, and global ones over a dataset.

We focus on the post-hoc global method of sensitivity analysis which describes how the model inputs influence the output (Dimopoulos, Bourret & Lek, 1995) by calculating partial derivatives of the model output with respect to its inputs. Since this method is model agnostic, it is also used as part of saliency maps in image classification with convolutional neural networks (Simonyan, Vedaldi & Zisserman, 2014), which show important parts of an image related to its classification. For more model agnostic post-hoc interpretability tools we refer the reader to Lundberg & Lee (2017) and Gevrey, Dimopoulos & Lek (2003).

The sensitivity analysis function in the prosper_nn package identifies the influences of the input variables on the output of PyTorch neural networks. The packages NeuralSens (Pizarroso, Portela & Muñoz, 2022) and NeuralNetTools (Beck, 2018) provide similar functionalities in R. In comparison to available implementations, our tool is not limited to feed forward networks but can handle recurrent structures as well.

Theory. A neural network can be seen as a function $f$ mapping an input vector $\mathbf{x}\in {\mathbb{R}}^n$ to an output $f(x)$. The partial derivative of the output with respect to an input variable $x_i$, $\frac{\mathrm{\partial }f(x)}{\mathrm{\partial }{x}_i}$, identifies the sensitivity of the output to this variable, i.e., how the output will change if the variable is altered. With a positive derivative the output will increase if the variable increases, and will decrease if the variable decreases. This relationship is reversed for negative derivatives. The absolute value of the derivation indicates the size of change to the output. Since these relations can be noisy, we propose to use the mean of an ensemble to have more reliable results (Mehringer et al., 2023).

If we calculate the sensitivities of all input variables for many observations, e.g., over time, we can classify the variables according to their influence on the model output:

• Constant: the derivations stay almost the same across the observations.

• Monotonic: the sign of the derivations remains the same across the observations.

• Non-monotonic: the sign of the derivations changes across the observations.

• Unrelated: the derivations are all close to zero.

These categories can be interpreted as follows: For a constant influence, there is a linear dependency between input and output. The direction of change depends on the sign. For monotonic results, the direction of change is fixed but the actual effect size can vary. If the sensitivity analysis shows a non-monotonic relation, the dependencies between input and output are more complex. For the unrelated sensitivity category, the output is not sensitive towards changes in the input variable. Hence, we can discard these variables for the purpose of feature selection.

Implementation. The sensitivity analysis function takes the trained ensemble and the data for which the sensitivity analysis is performed. Further, we set the output_neuron as a Tuple. From the tensor returned by the model, the tuple selects the node for which the analysis should be performed. This way, the function can process models with arbitrary output formats as long as they are a torch.Tensor and therefore also RNNs or CNNs.

sensitivity_analysis(model=ensemble, data=data,
    output_neuron=output_neuron)

The function returns the sensitivities of the output_neuron with respect to the input variables across the data observations graphically as a heatmap. Variables are displayed on the y-axis and observations on the x-axis. Negative relations are color coded in blue, positive ones in red. The opacity indicates the strength of the influence. Variables with constant influence reveal roughly the same shade of color over the observations. Monotone relations appear as either red or blue rows while non-monotone ones as rows of mixed color. Light colored rows identify variables with weak relations to the output.