Iván Cárdenas-Gallo

Raha Akhavan -Tabatabaei

Mauricio Sánchez-Silva

Emilio Bastidas-Arteaga

I Cárdenas
email: id.cardenas470@uniandes.edu.co

R Akhavan-Tabatabaei

A Markov regime-switching framework to forecast El Ni ño Southern Oscillation patterns

Keywords: ENSO, Markov-Switching, Autoregressive processes, Forecasting, Regimes 1 Introduction

The El Niño-Southern Oscillation (ENSO) is an ocean-atmosphere phenomenon involving sustained sea surface temperature fluctuations in the Pacific Ocean, causing disruptions in the behavior of the ocean and atmosphere. We develop a Markov switching autoregressive model to describe the Southern Oscillation Index (SOI), a variable that explains ENSO, using two autoregressive processes to describe the time evolution of SOI, each of which associated with a specific phase of ENSO. The switching between these two models is governed by a discrete time Markov chain (DTMC), with time-varying transition probabilities. Then, we extend the model using sinusoidal functions to forecast future values of SOI. The results can be used as a decision-making tool in the process of risk mitigation against weather and climate related disasters.

not completely predictable because of the complexity that results from the relationships between ocean currents, atmospheric circulation and winds in the Pacific. ENSO involves two cyclic phases; El Niño and La Niña. El Niño represents the warm phase of the ENSO. During this phase, the sea surface temperatures rise along the North-West coast of tropical South America. The impacts of El Niño episodes depend on the geographical location. For example, in Colombia El Niño results in long periods of droughts leading to problems in the agriculture sector, hydroelectric generation and water supply for vulnerable populations. El Niño is only one phase of the phenomenon. An episode of La Niña is followed most of the times by the El Niño episode. La Niña represents the cool phase of the ENSO and is associated with the cooling of ocean waters of the coast of Peru and Ecuador. Its impacts are completely opposite to those of El Niño. During La Niña conditions in Colombia, the amount of precipitations increases leading to floods, landslides and damages to civil infrastructure.

The importance of the phenomenon lies in its effect on the global climate. Although ENSO events are characterized by a fluctuation in ocean temperatures in the equatorial Pacific, they are also associated with changes in wind, pressure, and rainfall patterns. This set of climate changes cause a significant impact across the tropical Pacific region and areas away from it. Many studies have found a close relationship between Pacific sea surface temperature (SST) anomalies and precipitation (Li,2013). As mentioned above, during El Niño and La Niña the usual precipitation patterns can be greatly disrupted by either excessively wet or dry conditions (International Research Institute for Climate and Society, 2008) that may impact the global weather. Therefore, the ability to forecast El Niño and La Niña events is extremely important to define risk mitigation measures.

There are several indicators that provide information about the state of ENSO at a given time. One of them is the Southern Oscillation Index (SOI), which is calculated from the monthly fluctuations in the air pressure difference between Tahiti and Darwin. Sustained negative values of the SOI indicate El Niño episodes, while sustained positive values of the SOI indicate La Niña episodes. Another indicator of ENSO is the Multivariate ENSO Index (MEI), which is a monthly measure based on the six main observed variables over the tropical Pacific. MEI contains information of the sea level pressure, surface wind (East-West and North-South), sea surface temperature, surface air temperature and total cloudiness of the sky. In this case, sustained negative values of the MEI indicate La Niña episodes, and sustained positive values of the MEI indicate El Niño episodes. In Figure 1, we can observe the SOI and MEI from 1950 to 2012. Moreover, in Figure 1 it is possible to observe the two strong episodes of ENSO in the shaded areas: one El Nio episode from July-1993 to December-1993, and one La Nia episode from June-1973 to March-1976. [START_REF] Poveda | Hydro-climatic variability over the Andes of Colombia associated with ENSO: a review of climatic processes and their impact on one of the Earth's most important biodiversity hotspots[END_REF] reviewed the hydro-climatic variability of the Colombian Andes associated with El Niño-Southern Oscillation (ENSO) using records of rainfall, river discharges, soil moisture and a vegetation index (NDVI). Their conclusions point to define ENSO as the main forcing mechanism of interannual climate variability. Given the indexes that were described previously (SOI and MEI), they quantify the anomalies in the sea surface temperature in terms of their sign, timing and magnitude. These studies indicate that ENSO indexes become important and valuable tools to forecast several hydro-climatological variables in the region.

Over time, several models have been developed to forecast the behavior of the phenomenon. This effort has resulted in the atmospheric general circulation models (AGCMs). AGCMs of the atmosphere, ocean and sea ice were coupled and run asynchronously to produce credible simulations of the global climate (Meehl,1990). Given the implicit complexity to model the weather physically, this kind of approach becomes accessible only for national Most statistical models that have been developed for climatic variables are based on the Box-Jenkins methodology (BoxJenkins,1976). These models can represent the basic structure of the time series, but they fail to reproduce the nonstationary behavior that exists in climatic variables. A source of non stationarity is the existence of cyclic patterns that determine the behavior of the weather. For example, El Niño-Southern Oscillation (ENSO) is a cyclical process composed of a warm phase (El Niño) and a cool phase (La Niña).

Given the existence of cyclical processes, several models have been developed to include the fact of various changing regimes over time. For example, Ubilava et al. (2013) adopted a smooth transition autoregressive modeling framework to explain the sea surface temperature in order to forecast ENSO; they argued the advantage of nonlinear models to represent dynamic processes such as ENSO. In order to model the weather, the introduction of a variable that represents the weather type (regime) can be seen in [START_REF] Zucchini | Hidden Markov Models for Time Series: An introduction[END_REF] where Hidden Markov Models (HMM) are proposed for modeling the space-time evolution of daily rainfall.

HMMs have a variety of application, including the econometric models. [START_REF] Hamilton | A new approach to the economic analysis of nonstationary time series and the business cycle[END_REF] proposes econometric time series, which are composed of a Hidden Markov Model (HMM) and autoregressive models. In Markov Switching Autoregressive (MS-AR) models, several autoregressive models are used to describe the time evolution of the series and the switching between these different models is controlled by a Hidden Markov chain which represents the possible regimes. MS-AR models were extended by Andrew Filardo (Filardo,1994) who incorporated time-varying transition probabilities (TVTP) between regimes. Since its publication, several applications of the model have been implemented. Ailliot The main objective of this paper is to propose an MS-AR, in order to evaluate the pattern changes of ENSO along the time and to predict its behavior for the next years. This model contributes to the prediction of hydro-climatic variables, with significant practical implications for agriculture, hydropower generation, fluvial transport, natural hazards and disasters, and human health. The paper is organized as follows: Section 2 describes the proposed model analyzing its advantages and shortcomings. Section 3 presents the steps followed for the model implementation. Section 4 summarizes the results obtained showing the adequacy of the model to establish forecasts. Finally, several concluding remarks are given in the last section.

The Proposed -Markov Switching Autoregressive Model (MS-AR)

ENSO is a cyclical phenomenon in which it is possible to identify two regimes over time (La Niña and El Niño). The behavior patterns during each of the two phases are opposite and differentiable. On the other hand, the duration and the change between regimes is variable over time. Also, there is uncertainty about the length of episodes or phase transitions. This group of characteristics is properly evaluated by the Markov Switching Autoregressive Model (MSAR) proposed by [START_REF] Hamilton | A new approach to the economic analysis of nonstationary time series and the business cycle[END_REF] and improved by [START_REF] Filardo | Business-Cycle Phases and their transitional dynamics[END_REF]. In this section we present the model, its application and relevant considerations for its implementation.

Certain variables undergo episodes in which the behavior of the series seem to change quite dramatically according to a set of regimes. A MSAR model proposes a very tractable approach to model changes in regime. In order to develop a model that allows a given variable to follow a different time series process over time (y t ), consider a first-order autoregression in which the constant term c and the autoregressive coefficient f might be different for each regime s t proposed:

y t = c s t + f s t y t 1 + e t , (1)
where e t ⇠ i.i.d N(0, s 2 ). In order to model the change in regime over time we use a Discrete Time Markov Chain (DTMC). Then, we define s t as a random variable that can assume values according to a state-space {1, 2,...,N}. Also, the probability that s t+1 equals some particular value j depends on the past only through the most recent value s t :

P{s t+1 = j | s t = i, s t 1 = k....} = P{s t+1 = j | s t = i} = p i j (2)
The transition probability p i j gives the probability that state i will be followed by state j in a one-step transition. It is often convenient to collect the transition probabilities in an (N ⇥ N) matrix P known as the transition probability matrix.

In [START_REF] Hamilton | A new approach to the economic analysis of nonstationary time series and the business cycle[END_REF] assumes Fixed Transition Probabilities (FTP) for the DTMC. The assumption that the transition probabilities are time-invariant implies that the expected durations of phases do not vary over time. Due to the variability of the ENSO, this assumption might be too strong. El Niño occurs irregularly, with a recurrence ranging from two to eight years with variation in intensity. Each episode lasts about six to fifteen months. On the other hand, La Niña occurs less frequently than El Niño, and an episode lasts from nine months to three years. Therefore, we implemented a Markov-Switching estimation procedure that incorporates time-varying transition probabilities (TVTP) into the Discrete Time Markov Chain (Filardo,1994). There are advantages to using the MS-AR (TVTP) against MS-AR (FTP):

-The TVTP model allows the transition probabilities to decrease or increase after a change of regime. In this way, TVTP models have the flexibility to identify systematic variations in the transition probabilities after turning points. -The TVTP model captures more complex temporal persistence of the phases than an FTP model. Allowing the transition probabilities to change over time can represent long periods of a specific regime. -The TVTP model permits time varying expected durations of the phases. It represents a strong advantage over FTP, in which the expected durations of phases are constant along the time.

The TVTP Markov-Switching model for the values of SOI, y t is presented below:

y t = µ 1 + F(L)(y t 1 µ S t 1 ) + e t i f state 1 (Niño), (3)
y t = µ 2 + F(L)(y t 1 µ S t 1 ) + e t i f state 2 (Niña), (4)
where

F(L) = f 1 + f 2 L + ... + f r L r
where z t = {z t , z t 1 ,...,z t p } is referred to the values of MEI. In consequence, the monthly index of MEI becomes a support variable in the model. In the following section we present in detail the implementation of the proposed model.

Implementation

This section presents the implementation of the ENSO indexes, along with the description of our proposed methodology. First, it is necessary to identify the order of the autoregressive model (p) to capture the dynamic nature of the phenomenon. During this step, the autocorrelation function (ACF) and partial autocorrelation function (PACF) plots permit to identify the number of autoregressive terms that are needed, the estimation of the parameters including autoregressive coefficients and transition probabilities. Finally, we revise the model in order to determine whether the specified and estimated model is adequate. During the MS-AR model specification, the number of regimes and the order p of the autoregressive process were obtained. For our particular model, we defined two regimes trying to represent each phase of the ENSO. In this way, the state space of the Discrete Time Markov Chain is S t 2 {Niño(1), Niña(2)}. According to our experiments adding a third state associated with the neutral phase does not improve the results, while adding to the model complexity. The model can be observed in Figure 2. Moreover, the results obtained support the proposal of having only two states. On the other hand, we proposed an autoregressive process of order p = 2 based on the Bayesian Information Criterion (BIC). When choosing from several models, the one with the lowest BIC is preferred. Thus, the autoregressives processes are going to follow this structure:

y t = c 1 + f 1 1 y t 1 + f 1 2 y t 2 + e t i f state 1, (6)
y t = c 2 + f 2 1 y t 1 + f 2 2 y t 2 + e t i f state 2, ( 7
)
where f j i represents the autoregressive coefficient of order i for regime j, c j is the constant of the process for regime j 2 {Niño(1), Niña(2)}, and e t ⇠ i.i.dN(0, s 2 ) is a white-noise process.

For example, the dataset presented in Figure 1 that was used to build the model has a monthly record since 1950 to 2000 and was downloaded from the National Climatic Data Center (NCDC) in Colombia web site (National Climate Data Center, 2013). Given the model specification, the estimation of the autoregressive coefficients and the transition probabilities were developed in the computational software Eviews 8. The software uses the likelihood-based iterative algorithm proposed by [START_REF] Hamilton | A new approach to the economic analysis of nonstationary time series and the business cycle[END_REF], whose output is shown in Table 1. As it can be observed in Table 1, each one of the estimated coefficients is statistically significant (Prob. < 0.05). Consequently, the MS-AR(2) model is statistically significant to represent the data series of SOI. The standard errors and the p-values are estimated by the software EViews 8. The standard errors are computed using the generalized method of moments (GMM) formalized by [START_REF] Hansen | Large Sample Properties of Generalized Method of Moments Estimators[END_REF]. On the other hand, the p-values are obtained via the numerical method proposed by [START_REF] Hansen | The Likelihood Ratio Test Under Non-Standard Conditions: Testing the Markov Trend Model of GNP[END_REF]. Furthermore, the estimation of the transition probabilities required in the model was also carried out. Note that the transition probabilities change over time (TVTP); the results are presented in Figure 3. As it can be seen in Figure 3, the transition probabilities of being in the same state (p 11 , p 2,2 ) are approximated to one or zero for various periods. Each of these periods represent a historical episode of El Niño or la Niña. Thus, during the months in which an episode of El Niño occurred, the probability of remaining in regime 1 (p 11 ) is close to one. On the other hand, during these time intervals, the probability of remaining in regime 2 (p 22 ) is close to zero. For example, one historical event of El Niño took place during 1993. For this year, we can clearly observe the related behavior in Figure 3. As it might be expected, the same pattern can be observed during historical episodes of La Niña. A clear example could be the episode of 1975.

In Figure 4, historical episodes of El Niño are shaded in dark grey. During these time intervals, it can be seen that the probability of being in state 1 (El Niño) is approximately one. Likewise, historical episodes of La Niña are shaded in light grey. For these time periods, the probability of being in state 1 (El Niño) is close to zero. Based on these probabilities, we performed a validation process for the model. The duration of each of the episodes based on the model was identified and compared against the historical reports. The results obtained are presented in Table 2 From the transition probabilities, it is possible to infer the occurrence of El Niño and la Niña episodes. In this way, the MS-AR model allows to make inferences about when the episodes have occurred based on the observed behavior of the series. In Figure 4, the conditional probability of being in state 1 (El Niño) given the two previous observations of

y t (P[S t = 1 | y t , y t 1 , y t 2 ]
) is plotted. As we only have two states, the probability of being in state 2 (La Niña) is computed as the complement probability.

After the estimation of all the parameters of the MS-AR model, we verified its fitting to the data. In Figure 5 the actual values of SOI (1950 2000) are plotted against the values obtained by the proposed MS-AR(2) model. It can be observed that the proposed model captures the behavior of the index over time, highlighting the adequacy of the model. Also, it identifies the peaks and drastic changes in SOI along time. Likewise, it follows the trends The final stage of the implementation consists of verifying the fit of the model by analysing the residuals. The residuals are defined as the difference between the actual values of SOI and values obtained by the MS-AR(2) model proposed. Both set of values are plotted in Figure 5. The residuals should follow a white-noise process. Thus, it was assessed with statistical tests that the expected value of the residuals is zero. On the other hand, to verify that the residuals are not correlated, the corresponding autocorrelation function (ACF) was developed, and as a result of this we concluded that the residuals are not correlated for any lag. In the following section the implemented model is used to perform forecasts.

Forecasting SOI

One objective of this paper is to develop a stochastic model that allows us to predict future values of the SOI index, which is associated with ENSO. The MS-AR(2) model that was implemented in the previous section is useful to describe and analyze the data over time. However, an extension to the model is needed to forecast the SOI index. The main problem is the estimation of the transition probabilities of the Discrete Time Markov Chain. These are completely dependent on the support variable, for our particular application the MEI index. If we want to forecast future values of SOI, it is fundamental to estimate the transition probabilities for those months in the future. In this section we present a method to estimate future transition probabilities.

The proposed approach consists of representing the cyclic seasonality observed in the transition probabilities over time, with a sinusoidal function due to the cyclical behavior of the ENSO phenomenon. Thus, the time series (S t ) is represented as a periodic function with period s in which S t = S t s . The angular frequency w = 2p s indicates the angle in radians that is passed in a unit of time, taking into account that the complete cycle is 2p. To define the length of the period s, we used the expected durations of the El Niño and La Niña episodes. On average, El Niño epidodes last for 6 months, while La Niña epidodes have an average duration of two years. As a result, the length of the period s is established at 28 months. This implies that the transition probabilities p 11 t have a cyclical seasonality with s = 28 months. The time series was approximated with the following sinusoidal function.

p 11 t = µ + Rsin(wt + q ) + a t , (8)
where q is an initial offset angle that permits a more general adjustment to a data set. The constant R represents the amplitude of the cycle. Also, seasonal oscillations occur around the mean value of the series µ. Finally, we assume a random error a t with mean equal to zero, constant variance and normal distribution. To fit this model to the data series of p 11 , we modify Equation 8 using trigonometric identities:

p 11 t = µ + Rsin(wt)sin(q ) + Rcos(wt)cos(q ) + a t . (9)
Now defining A = Rsin(q ) and B = Rcos(q ), the expression concludes in:

p 11 t = b µ + b Asin(wt) + b Bcos(wt) + a t . (10)
The resulting model is linear with three unknown parameters which can be estimated by least squares. In this way, a sinusoidal function is fit to the monthly time series of the transition probabilities (p 11 ) which is shown in Figure 3.

An important contribution in the implementation of the MS-AR(2) model in Section 4 is related to the increase of persistence in time of the El Niño state. Over the years, the probability of being in state 1 (El Niño) seems to augment. To verify this behavior, we performed a 10 year moving average for the transition probabilities (p 11 and p 22 ). The result of this procedure is plotted in Figure 6. The increasing trend in the persistence of El Niño state is evident from Figure 6. This conclusion is supported by climatic studies that have been conducted to evaluate the effect of climate change on global phenomena. Results from several earlier global climate model studies indicated that as global temperatures rise due to increased greenhouse gases, the mean pacific climate will tend to resemble an El Niño -like state (Meehl,2000). This fact must be taken into account when estimating the transition probabilities. Consequently, we added the observed trend into the function described in Equation 10. Thus, the final function to estimate the transition probabilities is shown in Equation 14.

To evaluate the goodness of the predictions, the evaluation period is defined for the months between January 2000 and December 2005. Therefore, the transition probabilities were estimated for these months using the function given in Equation 15. The result of this procedure is shown in Figure 7.

Given the random nature of the Discrete Time Markov Chain, forecasts that are performed with the model are random variables. Monte Carlo simulation can be used to evaluate the accuracy of the model. For our particular application, we performed 1000 replications. Then, the monthly average values are computed and presented in Figure 8. The objective of the model is to capture patterns of behavior on the SOI index. The fact of forecasting exact values, peaks and drastic changes is not included in the scope of the model. Therefore, with the aim of properly comparing forecasts with historic values, the historic curve was smoothed with the technique of moving averages. The overall results of the forecast can be detailed in Figure 8. The error in the forecasts was measured using the mean absolute error (MAE = 0.6254) with respect to the smoothed curve. Moreover, given the random nature of the forecasts we computed 95% prediction intervals for each of the months in analysis. In Figure 9 the actual value of SOI and the lower and upper limits for each of the prediction intervals are plotted. For 58% of the months, the actual value of SOI was within the prediction interval.

The MS-AR(2) model is able to reproduce the behavior of the smoothed curve. When the index is negative, the average of the forecasts is also negative. This fact also happens in case of positive values. Therefore, it is possible to infer that the model is useful to predict the future behavior patterns of ENSO. The forecast for the next five years (December 2013 -December 2018) has been done and presented in Figure 10. As a result of the forecast done in November 2013, two episodes of El Niño are predicted to occur from September 2014 to June 2015 and July 2017 to January 2018. It should be noted that the predictions were performed with the average of the 1000 simulations. However, from the results it is possible to compute the probability distribution of the forecast in a histogram. For example, the histogram associated with the forecast of SOI for December 2014 is presented in Figure 11. In this way, each monthly forecast may be fit to a known distribution. This fact becomes more useful than obtaining a single value for each month because it is possible to observe the behavior of the forecast over time and calculate the probability that the index variable (SOI) exceeds a particular value. This information is useful to identify the probability of occurrence of extreme events.

Concluding remarks

This paper investigates the use of MS-AR model to describe time series associated with ENSO. Our methodology is based on a two regime Markov Switching Autoregressive model that allows separate estimation between El Niño and La Niña episodes. After implementing the model in Section 3 and extending it with a sinusoidal function in Section 4, forecasts were performed and compared with historical values. The obtained results indicate that the model fits well to the time series. On the other hand, the forecasts were able to capture the behavior pattern of the index SOI showing the virtues of the model. The proposed model is a very tractable approach to modeling time series with clear changes along the time. It is useful to describe variables to identify various regimes that the state of the weather fluctuates between. We consider this as the main advantage of our proposed method against the Box-Jenkins methodology that has been widely used in this field. Box-Jenkins methodology leads to a good description of the time series, but it does not capture certain non-linearities that occur due to the existence of changing regimes. Moreover, Box-Jenkins methodology does to perform forecasts in long-term. In the study of climatic variables, one source of nonlinearity in many meteorological time series is induced by the existence of weather cycles. Thus, the introduction of the cyclical phases of El Niño and La Niña allows to describe the ENSO patterns, properly.

The second advantage presented in the paper refers to the ability of the Markov Switching Model for defining and measuring cycles in a time series. With the model it is possible to draw probabilistic inference about whether and when they may occur based on the observed behavior of the series. This fact allows to discover regimes that are not directly observable and that could be used as explanatory factors of the variables in evaluation. Moreover, the interpretation of the model leads to open structure which allows more physical models. Such models consider the physical behavior of the phenomena to statistically explain the variables in analysis. On the other hand, the estimation of parameters by the method of maximum likelihood provides the foundation for forecasting future values of the series. As the main concluding remark, the Markov Switching Autoregressive Model is a methodology that can represent time series with changing regimes over time. In our particular application, the model explains the cyclical variability of the phenomenon properly, so predictions of future behavior patterns can be successfully performed. Finally, the most important contribution of the paper is that the results can be used as a decision-making tool in the process of risk mitigation against weather and climate related disasters. Future research directions include the implementation of this model with other climate variables and improvement of the model specifications to obtain better results.

Fig. 1

 1 Fig. 1 Evolution of the index SOI and MEI over time (1950-2012).

Fig. 2

 2 Fig. 2 Transition probability diagram of the DTMC including Time Varying Transition Probabilities.

Fig. 3

 3 Fig. 3 Evolution of the DTMC Transition Probabilities over time.

Fig. 4

 4 Fig. 4 Conditional probability of being in state 1 (El Niño) given the two previous observations of y t (P[S t = 1 | y t , y t 1 , y t 2 ]).

Fig. 5

 5 Fig. 5 Comparison between the MS-AR(2) model and the SOI index (1950-2000).

Fig. 6

 6 Fig. 6 DTMC transition probabilities 10 years moving average.

p

  11 t = b µ + b Asin(wt) + b Bcos(wt) + b tt + a t (11) p 11 t = 0.5043 0.2081sin(wt) 0.1202cos(wt) + 0.0016t + a t

Fig. 7

 7 Fig. 7 Estimated DTMC Transition Probabilities using Sinusoidal function.

Fig. 8

 8 Fig. 8 Comparison between the forecast obtained with the proposed model and the SOI index (2000-2005).

Fig. 9

 9 Fig. 9 Comparison between the actual value of SOI the 95% prediction intervals obtained with the proposed model.

Fig. 10

 10 Fig. 10 Forecast obtained with the proposed model for the SOI index (2014-2018).

Fig. 11

 11 Fig. 11 Histogram associated with the forecast of SOI for December 2014.

  Markov regime switching model. The hot and cold periods and their turning points are detected clearly (A hot period is related to positive market conditions). Fallahi et al. (2014) used Markov Switching models to identify and analyze the cycles in the unemployment rate and four different types of criminal activities in the US. Finally, Liu et al. (2006) developed a Markov regime Switching model with two regimes representing expansion and contraction to catch the cyclical behavior of the semiconductor business.

	et
	al. (2012) propose a Markov-Switching Autoregressive model to describe the behavior of
	wind time series obtaining a good fit to the data. Souza et al. (2010) used a Hidden Markov
	Model (HMM) to predict future crude oil price movements. Their results indicate that the
	proposed model might be a useful decision support tool. Walid et al. (2011) proposed a
	Markov Switching model to investigate the exchange rates for four emerging countries over
	the period 1994-2009. Yuan (2011) presents an exchange rate forecasting model which com-
	bines the multi-state Markov Switching model with smoothing techniques. Guo et al. (2010)
	focus on detecting hot and cold IPO (Initial Public Offering) cycles in the Chinese A-share
	market using a

  this way, a MS-AR process is a discrete-time process with two components {S

	t ,Y t }. t 2 {Niño(1), t denotes the monthly index of SOI with values in ( •, •) and S For our case, Y Niña(2)} represents the regime at month t. Then, the MS-AR process is characterized by the
	following two conditional independence assumptions:
	-The conditional distribution of S t given the values of {S t 0 } t 0 <t and {Y t 0 } t 0 <t only depends on the value of S t 1 . In this manner, we assume that the weather type S t is a first order
	Discrete Time Markov Chain.
	-The conditional distribution of Y t given the values of {S t 0 } t 0 <t and {Y t 0 } t 0 <t only depends on the value of S t and Y t 1 ,...,Y t p . In this manner, we assume that the value of SOI
	(Y

t ) is an autoregressive process of order p whose coefficients evolve in time according to the regime S t .

  [START_REF] Ailliot | Markov-Switching autoregressive models for wind time series[END_REF] is the lag polynomial, µ S t is the state-dependent mean, e t ⇠ i.i.d N(0, s 2 ) and S t 2 {Niño(1), Niña(2)}. The variation from the FTP model presented above appears in the transition probability matrix.

	P{S t+1 = s t+1 | S t = s t , z t } =		p 11 (z

t ) 1 p 11 (z t ) 1 p 22 (z t ) p 22 (z t ) ,

Table 1

 1

	Variable Coefficient	Std. Error	z-Statistic Prob.
			Regime 1 (El Ni ño)
	c 1	-1.052	0.192	-5.487	0.000
	f 1 1 f 1 2	0.428 0.199	0.065 0.068 Regime 2 (La Ni ña)	6.537 2.925	0.000 0.003
	c 2	0.640	0.165	3.889	0.000
	f 2 1 f 2 2	0.384 0.176	0.066 0.065	5.803 2.718	0.000 0.006

MS-AR(2) coefficients

Table 2

 2 . They indicate that the model adequately captures the duration of episodes of both regimes. Comparison between historical and the proposed model durations of ENSO episodes

	Episode	Historical Duration Model Duration
		Regime 1 (El Ni ño)
	March 1951 to February 1952	12	12
	March 1957 to November 1957	9	8
	September 1963 to March 1964	6	7
	March 1965 to November 1965	9	9
	June 1969 to February 1970	9	9
	March 1972 to December 1972	9	9
	June 1977 to November 1977	6	5
	April 1982 to February 1983	11	12
	May 1987 to January 1988	10	11
	July 1993 to December 1993	6	6
	March 1994 to December 1994	10	9
	April 1997 to March 1998	12	12
		Regime 2 (La Ni ña)
	January 1950 to February 1951	14	15
	April 1954 to January 1957	34	32
	April 1964 to October 1964	7	7
	June 1970 to March 1972	22	23
	June 1973 to Marc 1976	34	34
	April 1988 to July 1989	15	15


Joaquim Assunção
email: joaquim.assuncao@pucrs.br

Paulo Fernandes
email: paulo.fernandes@pucrs.br

Lucelene Lopes
email: lucelene.lopes@pucrs.br

Silvio Normey
email: silvio.gomez@pucrs.br

A DIMENSIONALITY REDUCTION PROCESS TO FORECAST EVENTS THROUGH STOCHASTIC MODELS

This paper describes a dimensionality reduction process to forecast time series events using stochastic models. As well as the KDD process defines a sequence of common steps to achieve useful information through data mining techniques, we propose a sequence of steps in order to estimate the probability of future events through stochastic modeling. Our process focus on reduce the dimensionality of data, thus reducing the effect of the common problems involved in stochastic modeling, such as the state space explosion and the large modeling efforts to create such models.

Introduction

Stochastic modeling can be a natural option to predict system behavior. This system can be a process that represents basic situations or even complex nature events. Either way, a good model can deliver interesting probabilities about system events for its past, or even its future.

However, modeling is a non trivial task that requires a specialist on the domain that have a good knowledge about the scenario to be modeled. Other problem frequently faced, with stochastic models, is the state space explosion that causes a limit of states, thus, a limit of data to be modeled.

The amount of data in the world, in our lives, seems to go on and on increasing [START_REF] Witten | Data Mining: Practical Machine Learning Tools and Techniques, Second Edition[END_REF], and as the volumes of data increase, it increases the difficulty of handling these data. One of the problems most frequently faced is the difficult to detect information among such large amounts of data. These informations are lost because it is hard for a human being to see the patterns and relations among such abundant data.

To cope with this problem, specialized data mining algorithms retrieve information that are practically invisible to our perception. Although the information is revealed, it is not yet knowledge, since knowledge must be obtained by rationally processing the extracted information. Aiming to formalize the steps commonly used to reach the knowledge through data mining, there are some well defined processes in the literature, as for example the knowledge discovery in databases (KDD) described by Han et al. [START_REF] Han | Data Mining, Second Edition: Concepts and Techniques, The Morgan Kaufmann Series in Data Management Systems[END_REF]. Although this process is well consolidated and often used as a guide for knowledge discovery, instead of describing the path to discover knowledge, it seems to have been created to cover the steps that usually come before and after data mining.

In the context of our paper we focus on a specific case of KDD, which is the forecast of future events in time series databases. Therefore, we propose the use of stochastic modeling techniques as steps of the proposed process, information about events is a common trade in stochastic models. Specifically, we simulate and analyze a formal state-based model of the system in order to retrieve probabilities for future events (and states) of the targeted system.

In a rough comparison. as well as KDD shows a set of steps commonly used to prepare data for data mining, here we describe a process that shows useful steps to prepare data for stochastic modeling. These steps aim to decrease the time spent to develop models, coping with state space problem and reducing the chance of human mistakes while manually developing complex models.

This paper is organized as follows: the next sections describe this paper background: Time series and Stochastic modeling formalisms; The fourth section puts this work in perspective with pre-existent techniques; The fifth section describes the proposed process; Section 6 illustrates the proposed process with some numerical result obtained applying it to pre-existent datasets; Finally, the conclusion summarizes our contribution and suggests future works.

Time Series

A popular method to represent and analyze data is time series (TS) [START_REF] Chakrabarti | Locally adaptive dimensionality reduction for indexing large time series databases[END_REF][START_REF] Lin | A symbolic representation of time series, with implications for streaming algorithms[END_REF][START_REF] Wang | Experimental comparison of representation methods and distance measures for time series data[END_REF]. Given the ease of use, modify, compare and analyze data; TS became very popular in many fields of sciences such as statistics, economy, biology, geography, seismology, engineering, communication, machine learning, etc. A TS is characterized as a collection of data spread in time. Formally, a TS S can be defined as follows:

S = [(p 1 ,t 1 ), (p 2 ,t 2 ), ..., (p k ,t k ), ...(p n ,t n )] (2.1)
where (t 1 < t 2 < ... < t k < ... < t n ) In this representation, p i is the data corresponding to the i-th observation, and t i is the time when p i was observed. In such way, a TS can be visualized in two axis, one with the data representing p i , and other with the time passage from t 1 to to t n . Despite of that, since the data associated to each observation can be be related to a complex information, TS are essentially high dimensional data [START_REF] Han | Data Mining, Second Edition: Concepts and Techniques, The Morgan Kaufmann Series in Data Management Systems[END_REF]. In fact, TS can be employed to describe very complex structures like fossil shapes, molecules, geological materials, etc. [START_REF] Ye | Time series shapelets: a new primitive for data mining[END_REF][START_REF] Keogh | Lb_keogh supports exact indexing of shapes under rotation invariance with arbitrary representations and distance measures[END_REF].

In order to enhance the ease of storage, to facilitate human analysis, and to be flexible to many domains, TS had became a global trending for representation of data. Consequently, the constant use of TS in different areas and the growth of machine learning, brought researchers to new challenges related to TS.

We can highlight such challenges according to the main goals of TS analysis: modeling and forecasting [START_REF] Han | Data Mining, Second Edition: Concepts and Techniques, The Morgan Kaufmann Series in Data Management Systems[END_REF]. In modeling, researchers try to improve efficiency by reducing the impact of the dimensionality factor and by improving measurements techniques. In forecasting, researchers try to predict future situations, e.g., weather forecast, traffic behavior, etc.

One of the big problems related to TS is how to improve modeling efficiency being TS generally a data streaming [START_REF] Wang | Experimental comparison of representation methods and distance measures for time series data[END_REF]. Aiming to solve this problem, in the last decade, many methods to represent TS have been created [START_REF] Chan | Efficient time series matching by wavelets[END_REF][START_REF] Chakrabarti | Locally adaptive dimensionality reduction for indexing large time series databases[END_REF][START_REF] Korn | Efficiently supporting ad hoc queries in large datasets of time sequences[END_REF].

Stochastic Modeling

Stochastic modeling is the art of representing the behavior of a system by describing its states and transitions. One of the most common among such formalisms are the Markov chains. Given an initial state, even a simple discretetime Markov chain can be used to forecast future states by achieving probabilities for each possible state for the next n steps [START_REF] Stewart | Probability, Markov Chains, Queues, and Simulation[END_REF]. These steps are given by the Chapman-Kolmogorov equation. Being P a probability matrix:

P (m) i j = ∑ all k P (l) ik P (m-l) k j f or 0 < 1 < m (3.2)
Being a classic formalism, Markov chains are easy to handle and use to perform quantitative analysis. Unfortunately large volumes of data tend to generate huge models that are hard to handle using an ordinary Markov chain (state space explosion). In such cases, the use of a structured formalism as Stochastic Automata Networks (SAN) [START_REF] Plateau | On the stochastic structure of parallelism and synchronization models for distributed algorithms[END_REF][START_REF] Brenner | The Need for and the Advantages of Generalized Tensor Algebra for Kronecker Structured Representations[END_REF] may reduce this problem impact.

SAN is a formalism that describes a complete system as a collection of subsystems that interact with each other in a Markovian behavior. The main advantage of SAN is that we can represent very large and complex system in a human understandable model that has an equivalent Markov chain model. Figure 1 exemplifies a SAN model with two automata and its equivalent Markov chain.

SAN formalism uses synchronized events and functional rates to build the system. Once a model is ready, with few automata it can be equivalent to a giant Markov chain. In addition, the tensor format of SAN makes possible the use of specialized algorithms [START_REF] Fernandes | Efficient descriptor-vector multiplication in Stochastic Automata Networks[END_REF][START_REF] Czekster | Split: a flexible and efficient algorithm to vectordescriptor product[END_REF][START_REF] Czekster | Efficient vector-descriptor product exploiting time-memory trade-offs[END_REF] which are more efficient than those used for solving large Markov chains [START_REF] Czekster | Split: a flexible and efficient algorithm to vectordescriptor product[END_REF][START_REF] Czekster | Efficient vector-descriptor product exploiting time-memory trade-offs[END_REF]. A complete SAN model is represented by a SAN code that can be executed [START_REF] Brenner | PEPS2007 -Stochastic Automata Networks Software Tool[END_REF] to generate a set of probabilities, including those which may be generated in a Markov chain.

Recently we successfully created a SAN model to represent and predict classes in geological events [START_REF] Assunção | A structured stochastic model for prediction of geological stratal stacking patterns[END_REF]; a different approach that returns results much similar to a data mining classification. Due to the processes of analysis and data collection, this model took much time to be created. However, even high dimensional data can be transformed into TS, and then, into a string that can be transformed both into a SAN model or into a Markov chain.

Dimensionality Reduction and Forecast

Both research areas, TS and stochastic modeling, have an interesting thing in common; they both have many applications that are useful for knowledge discovery. The use of TS for data representation and forecasting is presented in a myriad of works [START_REF] Chatfield | Time-series forecasting[END_REF][START_REF] Moghadam | A heuristic method for forecasting chaotic time series based on economic variables[END_REF][START_REF] Peralta | Shuffle design to improve time series forecasting accuracy[END_REF][START_REF] Linjessica | Neighborhood counting for financial time series forecasting[END_REF][START_REF] Box | Time series analysis: forecasting and control[END_REF]. Despite such abundant material, is also frequent to found the use of Markovian models to predict states in complex scenarios [START_REF] Engel | Can the markov switching model forecast exchange rates?[END_REF][START_REF] Yu | Short-term traffic flow forecasting based on markov chain model[END_REF][START_REF] Steffens | Real-time prediction in a stochastic domain via similarity-based data-mining[END_REF][START_REF] Brenner | Modelling Grid5000 point availability with SAN[END_REF][START_REF] Fernandes | Performance evaluation of software development teams: a practical case study[END_REF][START_REF] Assunção | A structured stochastic model for prediction of geological stratal stacking patterns[END_REF][START_REF] Fernandes | Analysis of exponential reliable production lines using kronecker descriptors[END_REF]. In this work, we use these both techniques trying to profit from the knowledge on TS modeling and stochastic prediction through Markovian models.

Using stochastic modeling techniques, we can retrieve accurate probabilities about a given system or event. However, the "state space explosion" makes it hard to scale models to real application size problems [START_REF] Valmari | The state explosion problem[END_REF][START_REF] Stewart | Probability, Markov Chains, Queues, and Simulation[END_REF]. In TS there is a similar phenomenon, the well-known "curse of dimensionality", which is the subject of many works in the area [START_REF] Chan | Efficient time series matching by wavelets[END_REF][START_REF] Chakrabarti | Locally adaptive dimensionality reduction for indexing large time series databases[END_REF][START_REF] Lin | A symbolic representation of time series, with implications for streaming algorithms[END_REF].

Actually, these two problems are different views of a same phenomenon, since the background limitation is the fact that it is easy to describe a large number of points of view (several dimensions), but the reality that needs to be handled is the Cartesian product of the data from all this points of view, i.e., we model focusing on the parts, but the model must be handle as a whole. Among the approaches to tackle this problem, we are interested in this paper in dimensionality reduction [START_REF] Lin | A symbolic representation of time series, with implications for streaming algorithms[END_REF] and structured modeling formalisms [START_REF] Brenner | The Need for and the Advantages of Generalized Tensor Algebra for Kronecker Structured Representations[END_REF].

Several techniques have been developed to represent TS in order to reduce the difficulty of processing high dimensional data. Here, we are interested in techniques that have two characteristic: 1) the representation must be symbolic; 2) the representation must be flexible in length. Among this techniques, for experiment purposes, we use Symbolic Aggregate approXimation (SAX) [START_REF] Lin | A symbolic representation of time series, with implications for streaming algorithms[END_REF].

SAX is a solution to reduce the dimensionality of a TS both in length (it reduces the number n of observations) and in data space (it aggregates the values of p i ). Also, it generates symbolic data that can be used to generate stochastic models on an automated way. The next section presents the process, and the basic phases to obtain the stochastic model, and then, the data forecasting through its transient and steady states analysis.

Proposed Process

Instead of work directly in mathematical models to slight improve some stochastic formalism, we focus on a process that reduces the data dimensionality before these data goes into the model. Instead of use a limited set of algorithms to forecasting in TS, we focus on adapt these data to work with advanced stochastic algorithms.

In order to achieve it, we carefully select the most useful techniques to represent and reduce the dimensionality of a given dataset, as well as, powerful methods that allow us to perform quantitative analysis under the generated data. Thus, this process aims to achieve a set of useful probabilities directly from a dataset using the minimum human effort and minimizing the effects of "state space explosion".

Figure 2 shows a simplified flow diagram of the proposed process. Each picture represented at the bottom part of this figure is a set of data; and each term written at the top part of this figure is a task of the process. Next, we describe each step of the process, which is basically composed by a task and its resultant data .

TS Representation

This step performs the dimensionality reduction of the dataset. One by one, the TS goes to SAX which performs a dimensionality reduction using Piecewise Aggregate Approximation [START_REF] Chakrabarti | Locally adaptive dimensionality reduction for indexing large time series databases[END_REF]. Then, SAX generates the symbolic data.

Coder

In big data, coding is the process of tagging terms using an identifier that is assigned to every synonymous term in the data [START_REF] Berman | Principles of Big Data: Preparing, Sharing, and Analyzing Complex Information[END_REF]. This step has a similar objective; although, instead of grouping terms, symbolic data are grouped into states. Thus, this step is the connection between TS output and the solution of the SM. In other words, Coder is the step responsible to get the symbolic data generated by TS representation and transform it into a SM. As seen in the Figure 3, the symbolic data is a collection of N symbols that represents N levels of a given TS; thus, represented by: S = [a, b, c, d, ...Ω] the symbolic data generated (D) can be formalized as: D i = {s ∈ S}. Once we have the symbolic data, it is possible to automatically create the SMs. Basically it is performed by getting the frequencies of changes for each class in a symbolic data.

Given three hypothetically TS, Q, S and T and its representing symbolic data w(Q), w(S) and w(T ), each one composed by classes from a to j. The coder calculates the frequency of changes for each possible combination, e.g.: a → b, a → c, ... , a → j b → a, b → c, ... , b → j ... , ... , ... , ... j → a, j → b, ... , j → j

Each frequency is turned into a percentage value that is used as a transition rate for the automata. Also, each set of symbolic data turns in one automaton. In this scenario, we will have a SAN composed by three automata created from w(Q), w(S), w(T). Having each one 10 states, e.g., a, b, . . . , j, our SAN model will be limited to 1000 reachable states, which is almost instantaneously to solve with current SAN algorithms.

Stochastic model solution

In this step, the code generated is used to create and solve a stochastic model. These stochastic models can be either, Markov chains, SAN or any other Markovian formalism. We develop a tool to make automatically create Markovian models from a symbolic data; furthermore, there is a work in progress to automatically generate SAN code [START_REF] Assunção | Unsupervised model generation for geological events[END_REF]. Once that we have a SAN model, SAN solvers like PEPS [START_REF] Brenner | PEPS2007 -Stochastic Automata Networks Software Tool[END_REF] and SAN Lite-Solver [START_REF] Sales | San lite-solver: a user-friendly software tool to solve san models[END_REF] can deliver steady state and transient solutions through specialized algorithms [START_REF] Fernandes | Efficient descriptor-vector multiplication in Stochastic Automata Networks[END_REF][START_REF] Czekster | Split: a flexible and efficient algorithm to vectordescriptor product[END_REF][START_REF] Czekster | Efficient vector-descriptor product exploiting time-memory trade-offs[END_REF].

Despite the high complexity of the stochastic model solution algorithms, it is possible to understand the way those solutions deliver their results considering the following simple example. Let us consider a stochastic model with states A,B and C. Given the probabilities: The main goal of the experiments is to measure model accuracy. Thus, aiming to avoid data bias, we chose different kinds of public datasets to test our process. A sample of four synthetic datasets were randomly picked from [START_REF] Keogh | The UCR Time Series Classification/Clustering Homepage[END_REF]: Coffe, Symbols, WordsSynonyms, and ItalyPwrDemand. A sample of four, real world, economy datasets were collected from Brazilian market activity stock [START_REF] Yahoo! | Brazilian market stock indices[END_REF]: petr3, bbas3, brfs3, and vale3. Paris Temp was collected from [START_REF] Datamarket! | The open portal to thousands of datasets from leading global providers[END_REF]; as a real world dataset, it represents the Paris monthly temperature from Le Bourge station. Finally, for a more controllable scenario, we made a synthetic dataset that describes a sound-wave reducing its amplitude over time: WaveDS.

For each dataset, we reserved 80% to training, and 20% to test the output. The probabilities is given by two basic analysis: transient states and steady state. We use 42 transient states to predict next values; we also use 10 symbols for each dataset and a variable value to the dimensionality reduction. This is specially useful because many datasets have a high density of values, in the TS Representation step, many near points in the X axis for a light change in the Y axis.

Considering the high impact of external factors, we need to create a controllable experiment aiming to check if the model is generating coherent results. For such task, we had to set an input that we already had an expected output. This dataset, called WaveDS, construction is depicted in Figure 4. Clearly, we can see that tending to infinity the value of the wave should be 0. The expected flow for this experiment is the sequence:

1. SAX reads the TS and returns the corresponding symbolic data, being "a" the class close to -20 and "j" the class close to 20.

2. The Coder reads the symbolic data, calculate the probabilities and creates the stochastic model.

3. The model must return a very high steady state percentage chance for the classes close to 0.

The three steps had worked as expected; the experiment returns a high percentage values for the steady state in classes "e" and "f" that represents values close to zero, 48.63% and 48.03%, respectively. This experiment shows that the expected sequence was performed and the process has worked as expected.

The same approach was performed for all other datasets. For the synthetic datasets the output match with the expected classes, although 2 in market datasets the expected class was the second and third with higher percentages. This was expected, because there are many fluctuations on the market datasets that are very difficult to predict using only the historical data.

From all datasets tested, we divided nearly 80% to run in our process and 20% to test the output. For the datasets Symbols and WordsSynonyms, we had no significant changes, i.e., the steady state for all classes was achieved in the very beginning, and the value became close to the original. In Paris Temp dataset, we found two high scores for non adjacent classes, i.e., classes that cannot be agglutinated into one because they are non contiguous. Although, for the other datasets our quantitative analysis was accurate with the test part of the dataset.

All experiments results are show in Table 1. In this table is indicated the name of the dataset, the probabilities of the more frequent classes, and a general indication if the experiment prediction was accurate (" " indicates an experiment with a correct prediction, i.e., the predicted class, " †" indicates otherwise).

Dataset

Predicted For dataset Paris Temp the problem seems related to the bimodal characteristic of the data, since two non adjacent classes (c and i) were the more frequent ones. Dataset vale3, however, seems to be a more complex case and it deserves a deeper analysis in the future.

Conclusion

We successfully created a process that drives us into a next step to forecast data with high dimensionally. Our process is capable of making stochastic predictions through Markovian models coping with the curse of dimensionality. Our tests shows that the new process created is useful to reduce the data dimensionality and perform forecasting through stochastic models. Also connecting the best techniques in the literature we were able to automate the steps that usually demands time and effort from specialists.

More than forecasting, this process provides a new way to perform data mining tasks in high dimensional data. Now we are able to use Makovian models to knowledge discovery with a large amount of data, due to the curse of dimensionality it was infeasible before. This brings us to new challenges, once that we have a large amount of algorithms to solve problems using stochastic models.

In the implementation front, the next step is to generate the SAN code through the Coder, thus we can save time performing multiples executions due to the SAN format and the large variety of the algorithms already implemented in the SAN tools [6] [27]. This certainly will bring us more performance to make more tests and use the model into new scenarios, such as big data.

Other possible application is the real time prediction for streaming data. A dimensionality reduction can be applied online generating SAN code to be solved in real time, thus predicting next states.

Figure 1 :

 1 Figure 1: SAN model example with 2 automata and its equivalent Markov Chain.

5. 1 Figure 2 :

 12 Figure 2: Process Phases

Figure 3 :

 3 Figure 3: Classes and symbolic data generated for Coffee dataset, an example of TS dataset.

2 6

 2 0.3 0.1 0.4 B. to 0.2 0.8 0.4 C. to 0.5 0.1 0.2 A transient solution can be derived applying Chapman-Kolmogorov Equation [29] we shall achieve the following probabilities to be in each state for the next step: A B C A. to 0.31 0.42 0.27 B. to 0.15 0.70 0.15 C. to 0.24 0.48 0.28 Finally the steady state solution is computed considering an infinite number of steps. In our simple example it results in the following values: A B C A. to 0.2 0.6 0.2 B. to 0.2 0.6 0.2 C. to 0.2 0.6 0.Experiments and Results

Figure 4 :

 4 Figure 4: Synthetic dataset, WaveDS.

Table 1 :

 1 Probabilities and classes found for each dataset (" " indicates accurate and " †" indicates inaccurate).

		probabilities Accuracy
	Coffee	j: 94.28%
	Symbols	f: 12.51%; g: 11.94%
	WordsSynonyms	i: 29.51%; j: 29.01%
	ItalyPowerDemand	i: 40.21%; j: 42.54%
	petr3	b: 18.68 %
	bbas3	b: 45.37%
	brfs3	b: 55.80%
	vale3	a: 37.30%	†
	Paris Temp	c: 21.55%; i: 25.08%	†
	WaveDS	e: 48.63%, f: 48.03%