{'fraction_non_alphanumeric': 0.03270596309690362, 'fraction_numerical': 0.00902779610907267, 'mean_word_length': 4.819354838709677, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Markovian RNN: An Adaptive Time Series Prediction Network with HMM-based Switching for Nonstationary Environments

We investigate nonlinear regression for nonstationary sequential data. In most real-life applications such as business domains including finance, retail, energy and economy, timeseries data exhibits nonstationarity due to the temporally varying dynamics of the underlying system. We introduce a novel recurrent neural network (RNN) architecture, which adaptively switches between internal regimes in a Markovian way to model the nonstationary nature of the given data. Our model, Markovian RNN employs a hidden Markov model (HMM) for regime transitions, where each regime controls hidden state transitions of the recurrent cell independently. We jointly optimize the whole network in an end-to-end fashion. We demonstrate the significant performance gains compared to vanilla RNN and conventional methods such as Markov Switching ARIMA through an extensive set of experiments with synthetic and real-life datasets. We also interpret the inferred parameters and regime belief values to analyze the underlying dynamics of the given sequences.


A. Preliminaries
We study nonlinear time series prediction with recurrent neural networks in nonstationary environments. In particular, we receive a sequential data and predict the next samples of the given sequence based on the knowledge about history, which includes the previous values of target variables and side information (exogenous variables). Time series prediction task is extensively studied for various applications in the machine learning [1]- [3], ensemble learning [4]- [6], signal processing [7]- [9], and online learning theory [4], [10] literatures. In most real-life scenarios such as finance and business applications, time-series data may not be an output of a single stochastic process since the environment can possess nonstationary behavior. In particular, the dynamics of the underlying system, which generates the given sequence can exhibit temporally varying statistics. Moreover, the behavior of the system can even be chaotic or adversarial [11], [12]. This  Therefore, successfully modeling the nonstationarity of the data carries importance while performing prediction.
Although linear models have been popular, partly since they have been integrated to most statistics and econometrics software packages, neural network-based methods are becoming widely preferred for the time series prediction task thanks to their ability to approximate highly nonlinear and complex functions [13]. In particular, deep neural networks (DNNs) with multiple layers have been successful in resolving overfitting and generalization related issues. Although their success in some fields such as computer vision has been demonstrated in numerous studies [14], this multi-layered structure is not suitable for sequential tasks since it cannot capture the temporal relations in time-series data properly [15]. To this end, recurrent neural networks (RNNs) are used in sequential tasks thanks to their ability to exploit temporal behavior. RNNs contain a temporal memory called hidden state to store the past information, which helps them to model time-series more successfully in several different sequential tasks [2], [16]- [19]. Hence, we consider nonlinear regression with RNN-based networks to perform time-series prediction.
To address the difficulties raised by nonstationarity, several methods mostly based on mixture of experts [20]- [23] are proposed due to their ability to represent nonstationary or piecewise sequential data. Mixture of experts models rely on the principle of divide and conquer, which aims to find the solution by partitioning the problem into smaller parts. These models usually consist of three main components: separate regressors called experts, a gate that separates the input space into regions, and a probabilistic model that combines the experts [6]. The fact that mixture of experts simplifies complex problems by allowing each expert to focus on specific parts of the problem with soft partitioning provides an important advantage while modeling nonstationary sequences [6]. However, these methods require training multiple experts separately, which disallows joint optimization and end-to-end training. In addition, these methods rely on enough diversity among experts such that each expert makes errors at different regions in the feature space. Otherwise, their performance compared to single expert models becomes negligibly better or even worse, if none of the experts can fit the data well enough [24].
In certain studies, simple linear regressors such as autoregressive integrated moving average (ARIMA) models are considered as experts, where each individual expert specializes in a very small part of the problem [20]. To perform gating operation between experts, several adaptive switching approaches such as Markovian switching and transition diagrams are widely preferred [21], [25]- [29]. For instance, Markov switching models are particularly popular techniques for the time series prediction task, particularly in econometrics literature [21], [29]- [33]. These approaches employ multiple linear regressors (or classifiers), where each regressor is responsible for characterizing the behavior of the time series in a different regime. The switching mechanism between these regimes is controlled by a first-order Markov chain. Markov switching model can capture more complex dynamic patterns and nonstationarity, especially when the assumption of the existence of different regimes with Markovian transitions hold. This model and its variants are applied in analyzing and forecasting business, economic and financial time series [21], [26]. For instance, these models are used to identify business cycles, which consist of several regimes such as expansion and recession states [31]. However, none of these methods consider nonlinear regression with recurrent neural networks, which limits their capability to capture complex temporal patterns. Our model, Markovian RNN can be interpreted as a temporally adaptive mixture of experts model, where the regime-specific hidden state transition weights inside the RNN cell are employed as experts and HMM-based Markovian switching performs the gating operation between regimes. In this way, Markovian RNN can detect different regimes and focus on each of them separately through learning separate weights, which enables our model to adapt nonstationarity while making predictions.
Although there exists a significant amount of prior work on the time series prediction task in nonstationary environments, we, for the first time in the literature, utilize the benefits of recurrent neural networks and HMM-based switching for nonlinear regression of nonstationary sequential data. In this study, we introduce a novel time series prediction network, Markovian RNN, which combines the advantages of recurrent neural networks and Markov switching. Our model employs a recurrent neural network with multiple hidden state transition weights, where each weight corresponds to a different regime. We control the transitions between these regimes with a hidden Markov model, which models the regimes as part of a Markov process. In this way, our model can capture the complex sequential patterns thanks to RNNs, and handle nonstationary with Markovian switching. We also optimize the whole network jointly at single stage. Our model can also be extended using different RNN structures such as long short-term memory (LSTM) [34] and gated rectified unit (GRU) [35] networks as remarked in Section III-B. Furthermore, we demonstrate the performance gains of Markovian RNN in the time series prediction task with respect to the conventional forecasting methods such as ARIMA with Markov switching [30] and neural network-based methods such as vanilla RNN [36]. We perform extensive experiments over both synthetic and real datasets. We also investigate the inferred regime beliefs and transitions, as well as analyzing forecasting error in terms of mean squared error.

B. Prior Art and Comparisons
A considerable amount of research has been conducted in machine learning, signal processing and econometrics literatures to perform time series prediction in nonstationary environments [1], [11], [21]. Although there are widely embraced linear methods such as autoregression (AR), moving average (MA), autoregressive integrated moving average (ARIMA) and their variants in the conventional time series prediction framework, these methods fail to capture complex temporal patterns [8], [37], since they cannot fully capture nonlinearity. There exists a wide range of nonlinear approaches to perform regression in the machine learning and signal processing literatures [1], [8], [19]. However, these earlier methods suffer from practical disadvantages related to computation and memory. Furthermore, they can perform badly due to stability and overfitting issues [7].
To overcome these limitations, neural network-based methods have been increasingly popular thanks to the developments in optimization and neural network literatures [37], [38]. Most of the recent studies adopt recurrent neural networks and its variants for sequential tasks. Certain studies have successfully applied recurrent neural networks for language, speech and text processing tasks [1], [39]. This approach has also been used for anomaly detection in temporal sequences and time series classification [40], [41]. In this study, we employ RNNs considering their power to capture complex nonlinear temporal patterns and generalization capability over unseen data.
There are several studies, which adopt mixture of experts based approaches for time series prediction. A mixture of ARMA experts is considered in [20] to obtain a universal approximator for prediction in stationary sequences. In the work, the authors interpret the mixture of experts as a form of neural network. Various studies have developed universal sequential decision algorithms by dividing sequences into independent segments and performing switching between linear regressors [20], [27], [28]. However, since these works utilize linear regressors as experts, they perform poorly in challenging scenarios, particularly when the temporal patterns in the given sequence is more complex. To remedy these issues, mixture of heterogenous experts is proposed in [22] for nonstationary sequences. Another study also employs nonlinear regressors as experts for stock price prediction task [23]. However, the nonlinear regressors employed in these studies have multilayered perceptron architectures without any temporally recurrent structure. This layered structure poorly performs in capturing time dependencies of the data due to its lack of temporal memory. Therefore, DNNs provide limited performance in processing temporal data and modeling time series [15]. Instead, we employ RNNs to handle the temporal patterns in time series data. In addition, we jointly optimize the whole network at single stage whereas mixture of experts models require separate training sessions for each expert.
Designing the gating model in mixture of experts approach is as crucial as choosing the expert models. For instance, authors employ randomized switching mechanism based on transition diagrams in [27], [28]. In another study, authors use a gating procedure based on a fuzzy inference system in [23]. However, most studies, especially in the business domain and finance literature, prefer Markovian switching based approaches since they express financial cycles more accurately [42]. Certain earlier variants of this approach such as Hamilton model [30], Kim-Nelson-Startz model [43] and Filardo model [32] have been commonly preferred and specifically designed for the tasks in business and finance applications [42]. However, these statistical methods are not flexible in terms of modeling, since they employ linear models with certain assumptions such as sequences with varying mean and variance. In another study, a mixture of linear experts approach with hidden Markov gating model is proposed and trained with expectation-maximization (EM) algorithm [21]. Similar approaches have been applied for anomaly detection tasks as well. For instance, in [44], authors develop an adaptive HMM with an anomaly state to detect price manipulations. Although Markovian switching-based methods are commonly used for sequential tasks in nonstationary environments, few of them consider nonlinear models, which are mostly simple multi-layer networks. In addition, they usually require multiple training sessions and cannot be optimized jointly. However, we introduce a jointly optimizable framework, which can utilize the benefits of nonlinear modeling capability of RNNs and adaptive Markovian switching with HMMs in an effective way.

C. Contributions
Our main contributions are as follows: 1) For the first time in the literature, we introduce a novel time series prediction model, Markovian RNN, based on recurrent neural networks and regime switching controlled by HMM. This approach enables us to combine the modeling capabilities of RNNs and adaptivity obtained by HMM-based switching to handle nonstationarity. 2) We use gradient descent based optimization to jointly learn the parameters of the proposed model. The proposed sequential learning algorithm for Markovian RNN enables us to train the whole network end-to-end at single stage with a clean pipeline. 3) Our model can prevent oscillations caused by frequent regime transitions by detecting and ignoring outlier data. The sensitivity of the introduced model can readily be tuned to detect regime switches depending on the requirements of desired applications, or with cross-validation. 4) Through an extensive set of experiments with synthetic and real-life datasets, we investigate the capability of Markovian RNN to handle nonstationary sequences with temporally varying statistics. We compare the prediction performance of the introduced model with respect to vanilla RNN and conventional methods with Markovian switching in terms of root mean squared error (RMSE) and mean absolute error (MAE). We also analyze the inferred regime beliefs and transitions to interpret our model.

D. Organization
The organization of the paper is as follows. We define the time series prediction task and describe the framework that uses recurrent neural networks and hidden Markov models in Section II. Then we provide the introduced model, switching mechanism, and sequential learning algorithm for Markovian RNN in Section III. In Section IV, we demonstrate the performance improvements of the introduced model over an extensive set of experiments with synthetic and real datasets and investigate the inferred regimes and transitions. Finally, we conclude the paper in Section V with several remarks.

II. PROBLEM DESCRIPTION
In this paper, all vectors are column vectors and denoted by boldface lower case letters. Matrices are denoted by boldface upper case letters. x T and X T are the corresponding transposes of x and X. x is the 2 -norm of x. and denotes the Hadamard product and division, i.e., element-wise multiplication and division, operations, respectively. |X| is the determinant of X. For any vector x, x i is the i th element of the vector. x ij is the element that belongs to X at the i th row and the j th column. sum(·) is the operation that sums the elements of a given vector or matrix. δ ij is the Kronecker delta, which is equal to one if i = j and zero otherwise. Enotation is used to express very large or small values such that mEn = m × 10 n .
We study nonlinear prediction of nonstationary time series. We observe a vector sequence where T is the length of the given sequence, and x t ∈ R nx is the input vector for the t th time step. Input vector can contain target variables (endogenous variables) as well as side information (exogenous variables). The target output signal corresponding to x 1:T is given by y 1:T = {y t } T t=1 , where y t ∈ R ny is the desired output vector at the t th time step. Our goal is to estimate y t using the inputs until the t th time step bŷ where f is an adaptive nonlinear function parameterized with θ. After observing the target value y t , we suffer the loss (y t ,ŷ t ), and optimize the network with respect to this loss. We evaluate the performance of the network by the mean squared error obtained over the sequence with where MSE (y t ,ŷ t ) = e T t e t , and e t y t −ŷ t is the error vector at the t th time step. Other extensions are also possible such as mean absolute error (MAE) as remarked in Section III-B.

A. Recurrent Neural Networks
We particularly study time series prediction with RNNs [45]. For this task, we use the following form: where Here, h t ∈ R n h is the hidden state vector at time step t. W hh ∈ R n h ×n h , W xh ∈ R nx×n h and W hy ∈ R n h ×ny are the weight matrices. We use θ h = {W hh , W xh } and θ y = {W hy } to denote the state transition and state-to-observation parameters respectively.
We note that the introduced framework can be applied for any stateful neural network structure. Hence, it can be extended to various RNN-based networks such as Gated Recurrent Units (GRU) [35] and Long Short-term Memory Units (LSTM) [34]. We provide the equations for possible extensions in Remark 1 in Section III-A. Here, we consider RNNs due to its simplicity and practicality in real-life applications. We also do not state bias terms explicitly since they can be augmented into input vector such that x t ← [x t ; 1].

B. Hidden Markov Models
We utilize HMMs to model the switching mechanism of RNNs, as will be described in Section III. HMM is a statistical model, which consists of a discrete-time discrete-state Markov chain with unobservable hidden states k t ∈ {1, ..., K}, where K is the number of states. The joint distribution has the following form: The observation model (emission model) is sometimes modeled as a Gaussian such that p(y t |k t = k; θ) = N (y t |µ t , Σ t ).
The state posterior p(k t |y 1:T ) is also called the filtered belief state and can be estimated recursively by the forward algorithm [46] by ..] T as the K-dimensional likelihood vector respectively. Then (5) can be expressed as The filtered belief state vector can be obtained after normalizing the expression in (6) through dividing by the sum of values. We note that we call HMM states as regimes from now on to prevent terminological ambiguity with the hidden states of RNN. In the following section, we introduce the Markovian RNN architecture with HMM-based switching between regimes. We also provide the equations and the sequential learning algorithm of our framework.

III. A NOVEL RNN STRUCTURE
In this section, we introduce our novel contributions for sequential learning with RNNs. We provide the structure of the Markovian RNN, by describing the modified network with multiple internal regimes in Section III-A and HMMbased switching mechanism in III-B. We present the sequential learning algorithm for the introduced framework in Section III-C. The detailed schematic of the overall structure of our model is given in Fig. 1.

A. RNNs with Multiple Internal Regimes
Here, we describe the introduced Markovain RNN structure with multiple regimes, where each regime controls state transition independently. To this end, we modify the conventional form given in (2) and (3) as where k ∈ {1, ..., K} is the regime index, and K is the number of regimes. We also illustrate the modified RNN cell with multiple regimes in the left hand side of Fig. 1.
Here, the hidden state vector is independently propagated to the next time step at each node with different weights θ (k) h . We highlight that the state-to-observation model is same for all regimes. However, the resulting predictionsŷ (k) t are still different for each regime because of different hidden states h (k) t obtained for the t th time step. We obtain the final estimate of the hidden state at time step t by the weighted average of hidden states of each regime as where w t,k is the weight for the k th regime. Finally, we estimate the output using (3). Here, the weights w t,k are determined by the switching mechanism described in Section III-B. The number of states, K, is considered as a hyperparameter and can be selected using cross-validation.
Remark 1. Our model can also be extended with different RNN structures such as long short-term memory (LSTM) [34] and gated rectified unit (GRU) [35]. For instance, for LSTM, all gating operations and state updates can be performed for each internal regime with the following equations: are the forget, input and output gates, andc (k) t is the candidate cell state at time step t for the k th regime such that where σ and f h are nonlinear element-wise activation functions. For the final estimates of the hidden state, we can apply (9). For the cell state, the same form is applicable as well: The final output estimateŷ t can be calculated with (3).

B. HMM-based Switching Mechanism
We employ an HMM to control the switching mechanism between internal regimes. In particular, we perform soft switching, where the weight given in (9) are represented using the belief values of the HMM as follows: where α t−1,k w t,k denote the belief for the k th regime. To perform belief update as given in (6), we need to calculate the likelihood values of φ t for the t th time step after observing y t . To this end, for mean squared error loss, we consider the error model with Gaussian distribution such that where e (k) t is the error vector and R (k) t−1 is the error covariance matrix for the k th regime, which stores the errors of the corresponding regime until the t th time step, excluding the last step. Then we compute the likelihood by Once we obtain the likelihoods, we update the regime belief vector using (6) asα where we calculate φ t = [..., p(y t |k t = k), ...] T with (13). We finally update the error covariance matrix using exponential smoothing by where β ∈ [0, 1) controls the smoothing effect, which can be selected using cross validation. For instance, β = 0.95 would result in high sensitivity to errors, which can cause outlier data to bring frequent oscillations between regimes, whereas very small values for β might prevent the system to capture fast switches. The second part of the schematic in Fig. 1 illustrates the operations of HMM-based switching module.

Remark 2.
Our frameworks can also be used with different loss functions such as mean absolute error (MAE) loss. In this case, we can model the distribution of the error vector e t with the multivariate Laplacian distribution such that e t ∼ L(0, Σ t ), where Σ t is the error covariance matrix at the t th time step. The likelihood computation of regimes given in (13) can be modified for the multivariate Laplacian distribution as t , v = 1 − n y /2 and K v is the modified Bessel function of the second kind [47]. For onedimensional case (n y = 1), considering scalars instead of vectors, the likelihood equation reduces to where e k t and r k t are the error value and error variance at the t th time step for the k th regime.
Remark 3. HMM-based switching inherently prevents instability due to the frequent oscillations between regimes or possible saturations at well-fitted regimes. One might argue that certain regimes that have been explored heavily during the early stages of the training would dominate other regimes and cause the system to degenerate into quite a few number of regimes. However, since the error covariance matrix penalizes errors for well-fit regimes more harshly than the regimes that are not explored yet, the model will tend to switch to other regimes if the predictions made by the dominating regimes start to produce high errors. Here, the choice of the smoothing parameter β can be interpreted as an adjuster of the tolerance for the errors made in different regimes. As β increases, the switching mechanism will have greater sensitivity to errors, which can cause instability and high deviations in the regime belief vector. Likewise, as beta approaches towards zero, the system will not be able to capture switches due to the smoothing effect. This can eventually lead to saturations at well-fitted regimes. Thus, the choice of β directly affects the behavior of our model and we can readily tune it depending on the needs of the specific application or select it with crossvalidation. We further discuss and illustrate the effect of this parameter in Section IV-C.

C. Sequential Learning Algorithm for Markovian RNN
In this section, we describe the learning algorithm of the introduced framework. During training, at each time step t, our model predicts the hidden state h t and outputŷ t . We receive the loss given in (1) after observing the target output y t . We denote the set of weights of our model as We use the sequential gradient descent algorithm [48] to jointly optimize the weights during the training.
In Algorithm 1, we present the sequential learning algorithm for the introduced Markovian RNN. First, we initialize the model weights θ, hidden state h 1 , regime belief vector α 1 and error covariance matrices {R (k) 1 } K k=1 . For a given sequence with temporal length T , after receiving input x t at each time step t, we compute hidden states for each internal regime using (7). Then, we predict the output with these hidden states for each regime using (8). After forward-pass of each internal regime, we generate h t and output predictionŷ t using (10) and (3). After receiving the target output y t , we compute the loss using (1) and update the model weights through the backpropagation of the derivatives. We use the Truncated Backpropagation Through Time (TBPTT) algorithm [49], and provide the derivatives of model weights in Appendix A. To satisfy the requirement that each row of Ψ should sum to one, we scale the values row-wise using softmax function such that . Finally, we update the regime belief vector and error covariance matrices by (13)- (15).

26:
HMM Based Switching: return θ best 49: end if 50: end for 51: return θ best part, we simulate three synthetic sequences with two regimes and analyze the inference and switching capacity of our model under different scenarios. In the second set of experiments, we demonstrate the performance enhancements obtained by Markovian RNN in three real-life financial datasets and compare our results with the results of other methods. Also, we investigate the inferred regimes for these datasets by interpreting the temporal evolution of belief state vector and switching behavior of Markovian RNN.

A. Synthetic Dataset Experiments
In order to analyze the capability of Markovian RNN to detect different regimes, and to investigate the switching behavior between these regimes, we conduct initial experiments on synthetic data. We first describe the simulation setups for synthetic data generation and then, present the results obtained by all methods on the synthetic datasets. We also investigate the effect of error covariance smoothing parameter (β) introduced in (15) and mentioned in Remark 3.
1) Simulation Setups: In the synthetic data experiments, our goal is to predict the output y t given the input data x 1:t such that x t ∈ R is a scalar. The output data is given by y t = x t+1 . Here, the goal of these experiments is to conceptually show the effectiveness of our algorithm. In order to demonstrate the learning behavior of our algorithm with different patterns and switching scenarios, simulated sequences should have various regimes, where each regime possesses different temporal statistics. To this end, we conduct three experiments in which we simulate autoregressive processes with deterministic and Markovian switching, and a sinusoidal with Markovian switching.
a) Autoregressive Process with Deterministic Switching: In the first synthetic dataset experiment, we aim to generate a sequence with sharp transitions and obvious distinctions between regimes. To this end, we generate an autoregressive (AR) process with deterministic switching, which is given with the following set of equations: where x t ∈ R is the value of the time series at the t th time step, and ∼ N (0, 0.01) is the process noise. Here, (16) describes an AR process with two equal-duration (500) regimes, in which the system osciallates between. The first regime describes a random walk process, whereas the second regime gradually drifts towards white noise. The simulated system is deterministic in terms of switching mechanism between regimes since it strictly depends on the time step. Fig. 2a demonstrates the time series generated by this setup. b) Autoregressive Process with Markovian Switching: In this setup, we consider Markovian switching instead of deterministic switching. Here, the transition between regimes has Markovian property, therefore the regime of next time step only depends on the current regime. We consider third order AR processes with the coefficients of {0.95, 0.5, −0.5} and {0.95, −0.5, 0.5} for each regime respectively, and ∼ N (0, 0.01). For the transition matrix, we consider Ψ =  c) Sinusoidal with Markovian Switching: In this experiment, we generate a noisy sinusoidal signal with two regimes, where every regime represents a different frequency. Here, the simulated signal has two regimes with the magnitude of 0.5 and periods of 50 and 200 for the generated sinusoidals. The whole sequence consists of 5000 time steps and Markovian switching is controlled by the transition matrix Ψ = 0.99 0.01 0.01 0.99 . We also scale the magnitude to half and add Gaussian noise to the generated signal ( ∼ N (0, 0.0025)).
2) Synthetic Dataset Performance: Here, we present the training procedure and the results of the methods in terms of the root mean squared error (RMSE) and mean absolute error (MAE). In these experiments, each synthetic time series has 5000 time steps of temporal length. We split the data  into three splits for training (60%), validation (20%), and test (20%) respectively. We perform training on the training set and choose the best configuration based on the performance in the validation set. Then, we compare the test results of the best configuration for each method. In Table I, we provide the test RMSE and MAE obtained by each method on the synthetic datasets. In all setups, our model performs significantly better than other methods. Regardless of the switching mechanism and process dynamics, our model brings considerable improvements in terms of prediction RMSE and MAE. The performance Kim, Nelson and Startz model is not competitive with other methods, since it relies on switching variance between regimes. TVTP also has the form of MS-ARIMA, but it also models the transition probabilities as temporally varying values. In our setups, the transition probabilities are fixed, hence this method does not bring any improvement over MS-ARIMA. MS-ARIMA works significantly more accurately than standard ARIMA as expected. Likewise, Markovian RNN enjoys the benefits of adaptive HMM-based switching, which improves the predictions compared to the predictions of vanilla RNN.

B. Real Dataset Experiments
In this section, we provide the performance of our model and other methods in two financial datasets. We consider USD/EUR, USD/GBP and USD/TRY currency ratio data [50] from January 1, 2010 to January 24, 2020 with hourly resolution for these experiments. We calculate the mean and variance of prices for every day, and work in daily resolution. The goal is to predict the average currency rate for the t th day (y t ) given the input data x 1:t . In this case, the input data is a twodimensional vector such that x t ∈ R 2 since we consider x t to include the mean and standard deviation of the currency rate for the t th day. We split the data into three splits for training (60%), validation (20%), and test (20%) respectively.
For We take the first-order difference of the data to decrease the trend effect. Before testing, we calibrate the predictions by fitting a linear regressor on validation set, which minimizes the residual sum of least squares between validation predictions and validation labels. We also perform early stopping based on validation error such that we stop the training if the loss does not decrease for 20 consecutive epochs or the number of epochs reaches to 200. We obtained the best validation results for Markovian RNN in these experiments using the following parameters respectively: In Table II, we provide the test RMSE and MAE values obtained by each method on these currency datasets. All methods perform worst in USD/TRY due to the high oscillations in the sequence. In all cases, Markovian RNN performs significantly better than other methods. Our model obtains the lowest RMSE and MAE values thanks to the efficient combination of nonlinear modeling capability of RNNs and adaptive switching controlled by HMM. Although the regimes that affect the sequential dynamics cannot be observed, our model is able to detect these regions and switch between them to adjust the predictions for nonstationarity. We further analyze the switching behavior of Markovian RNN in Section IV-C and investigate the inferred regime belief vectors.   Fig. 3a, background colors represent the real regime value, where red color is used for the first regime and blue color is used for the second regime. Our model can properly distinguish between the two regimes, thus the resulting belief state vector values are consistent with the real regimes. In Fig. 3b, since the experiment is performed on a real-life dataset, in which the real regime values are not observable, consistency analysis is not possible. However, we still observe that our model switches between regimes in a stable way without any saturation. .

C. Regime Switching Behavior of Markovian RNN
To investigate how our model detects the regimes and decides to switch from one to another, we illustrate the belief state vector through time in Fig. 3a and 3b for AR process sequence with deterministic switching and USD/EUR dataset, respectively. For instance, in Fig. 3a, we observe that our model can determine different regimes correctly and adjust the belief vector accordingly. Here, background colors are used to indicate the real regime. During the analysis of reallife dataset experiments, checking the consistency between belief state vector and real regimes is not possible, since they are not observable. In Fig. 3b, we observe that there are certain periods in which the second regime dominates the predictions. Considering the switching behavior in Fig. 3b, we can interpret that our model does not suffer from rapid oscillations between regimes or does not saturate at a regime. The evolution of the belief state vector is stable, but still responsive to nonstaionarity.
In addition, we zoom in to the predictions of Markovian RNN around switching regions in Fig. 4a, and show that our model can successfully adapt to nonstationarity. As we mention in Remark 3, the error covariance smoothing parameter (β) in (15) can be tuned with cross-validation to adjust the switching behavior of Markovian RNN. For instance, lower values of this parameter (β = 0.5) bring more toleration for prediction errors and may provide robustness against outliers with the expense of lagged transitions between regimes as shown in Fig. 4a. On the contrary, higher values (β = 0.9) provide faster transitions but may cause oscillating predictions or lead to instability.

V. CONCLUSION
We study nonlinear regression for time series prediction in nonstationary environments. We introduce a novel time series prediction network, Markovian RNN, which is an RNN with multiple internal regimes and, HMM-based switching.
Each internal regime controls the hidden state transitions with different weights. We employ an HMM to control the switching mechanism between the internal regimes, and jointly optimize the whole network in an end-to-end fashion. By combining the nonlinear representation capability of RNNs and the adaptivity obtained thanks to HMM-based switching, our model can capture nonlinear temporal patterns in highly nonstationary environments, in which the underlying system that generates the time series has temporally varying dynamics.
Through an extensive set of synthetic and real-life dataset experiments, we demonstrate the performance gains compared to the conventional methods such as vanilla RNN [36], MS-ARIMA [30] and Filardo model [32], which are commonly preferred in business, economy and finance applications. We show that the introduced model performs significantly better than other methods in terms of prediction RMSE and MAE thanks to the joint optimization and the efficient combination of nonlinear regression with RNNs, and HMM-based regime switching. We show that Markovian RNN can properly determine the regimes and switch between them to make more accurate predictions. We also analyze the effect of the error covariance smoothing parameter on the switching behavior of our model. As the experimental results and our analysis indicate, our model can capture nonlinear temporal patterns while successfully adapting nonstationarity without any instability or saturation issues.

APPENDIX A
In this part, we provide the derivatives of the model weights RNN. The equations of the basic derivatives are as follows: ∂α t,k ∂α t,k = δ kk sum(α t ) −α t,k sum(α t ) 2 , where z xh x t . We use (19) and (20) to obtain the following: where ∂ht ∂ht−1 = K k=1 α t−1,k W (k) hh diag(f h (z (k) t )). Then, we can use (17)- (20) and (23) to obtain ∂ t ∂W (k) xh and ∂ t ∂W (k) hh as follows: where t ), and τ is the truncation length. Here, d is a vector such that d i = δ ii .