Bayesian Analysis of Multiplicative Seasonal Threshold Autoregressive Processes

Seasonal fluctuations are often found in many time series. In addition, non-linearity and the relationship with other time series are prominent behaviors of several, of such series. In this paper, we consider the modeling of multiplicative seasonal threshold autoregressive processes with exogenous input (TSARX), which explicitly and simultaneously incorporate multiplicative seasonality and threshold nonlinearity. Seasonality is modeled to be stochastic and regime dependent. The proposed model is a special case of a threshold autoregressive process with exogenous input (TARX). We develop a procedure based on Bayesian methods to identify the model, estimate parameters, validate the model and calculate forecasts. In the identification stage of the model, we present a statistical test of regime dependent multiplicative seasonality. The proposed methodology is illustrated with a simulated example and applied to economic empirical data.


Introduction
The threshold autoregressive models (TAR) proposed by Tong & Chen (1978) and Tong & Lim (1980) have been widely applied in several areas, including economics, finance, hydrology, meteorology and biology. The TAR models assume that the values of a stochastic process {Z t } (threshold process) determines both the values of the stochastic process {X t } (interest process) and its dynamics. When the threshold process is the same process of interest but lagged, the model is known as SETAR (Self-Exciting TAR), Tong (1990). The TAR models can be easily extended, including some exogenous time series and called TARX. The use of Bayesian methods applied to the TAR and SETAR models, can be found in Chen & Lee (1995), So & Chen (2003), Nieto (2005), Nieto (2008), Nieto, Zhang & Li (2013), Zhang & Nieto (2015) and Calderón & Nieto (2017). While to the TARX models, can be found in Chen (1998), So, Chen & Liu (2006), , Chen, Gerlach & Lin (2010), among other works. Aditionally, Hansen (2011) and Chen, Liu & So (2011) give a review of the theory and applications of the TAR and TARX models in economics and finance, respectively. TAR models may be useful for nonseasonal or seasonally adjusted data that display nonlinearity by means of thresholds. However, diverse time series present some seasonal patterns. In addition, the nonlinearity of thresholds and the relationship with other time series are typical characteristics of many such time series in practice. Tong (2015) gives a compendium of what has been developed with respect to the TAR models, and possible issues to be investigated, including the study of seasonality. In many cases, it is convenient to study the seasonal patterns present in time series, since they can carry relevant information regarding the dynamic behavior of the series.
Thus, using a frequentist approach and data without seasonal adjustment, Franses, de Bruin & van Dijk (2000) and Franses & Van Dijk (2005) examine the performance of various models of linear and smooth transition autoregressive (STAR) time series, considering deterministic seasonality of series of quarterly industrial production data. Crespo (2001) considers seasonality as deterministic and dependent of the regime in the SETAR models and uses his model in time series of quarterly unemployment rates. De Gooijer & Vidiella-i Anguera (2003) define a seasonal SETAR model, where each regime has the form of a multiplicative seasonal autoregressive model. They carry out the study of time series of monthly inflation rates in several countries, finding more realistic dynamic characteristics than those explained by the usual SETAR or linear seasonal time series models. In this paper, we propose a statistical model called the multiplicative seasonal threshold autoregressive model with exogenous input (TSARX), which explicitly and simultaneously incorporates multiplicative seasonality by regimes that arises from periodic autocorrelation structures and threshold nonlinearity. The proposed model arises from imposing certain restrictions on the parameters of a TARX model. Bayesian methods are used to determine the presence of multiplicative seasonality by regimes and to obtain joint estimates of parameters, generalized residuals and forecasts. The developed methodology is illustrated with a simulated TSARX model and with an empirical application in the Colombian economy. The variables to be used are the growth rates of total monthly unemployment and a monthly coincident index that traces the Colombian economy (ISE for its Spanish acronym of índice de seguimiento económico).
The paper is organized as follows: In Section 2, the TSARX model and the proposed methodology are presented. Section 3 shows a simulated example. The real-data application is presented in Section 4 and, finally, Section 5 concludes.

Model Specification
Let {X t } and {Z t } be stochastic processes related by means of the equation if r j−1 < Z t−d ≤ r j , for all t ∈ Z (Z is the set of integers numbers) and some j = 1, . . . , l. The model assumes l regimes defined by the threshold values r j that satisfy −∞ = r 0 < r 1 < · · · < r l−1 < r l = ∞ and the threshold variable lagged d periods Z t−d . In addition, the threshold variable is considered an exogenous input, in the sense that there is no feedback from {X t } towards {Z t }. The parameters of the model can be divided into two groups (Nieto 2005): • Structural parameters: the number of regimes l, the l − 1 thresholds r = (r 1 , . . . , r l−1 ) ′ , the delay parameter d, the autoregressive orders of the l regimes k 1 , . . . , k l , K 1 , . . . , K l , q 1 , . . . , q l and seasonal period s.
We assume that the seasonal period of model s is known, where s = 4 is usually taken for quarterly time series and s = 12 for monthly time series. Moreover, model (1) could be used for time series that exhibit another type of seasonal patterns, as, for example, either daily (hour data) or weekly (daily observations). We consider k j < s for all j(j = 1, . . . , l).
Also, {ϵ t } is a zero mean Gaussian white noise process with variance 1 (GW N (0, 1)) such that E(X w ϵ t ) = 0 for w < t, Z w and the set {X t , X t−1 , . . .} are mutually independent for all w > t and {Z t } and {ϵ t } are mutually independent. To describe the stochastic behavior of {Z t }, we additionally assume that {Z t } is a homogeneous Markov chain of order m ≥ 1, with initial density function f (·) and kernel density function f m (· | ·), in the Lebesgue -measure sense. Additionally, the Markov chain converges weakly to a distribution F . For more details about Markov chains with general state space, see Meyn & Tweedie (2009). Equation (1) describes a dynamic system without feedback, with observable input {Z t } and observable output {X t }, which is presented in Tong (1990) but without considering the modeling of multiplicative seasonality by regimes.
The symbol TSARX(l; k 1 , . . . , k l ; K 1 , . . . , K l ; q 1 , . . . , q l ) s , is used to denote model (1) and it is said that {X t } is a threshold multiplicative seasonal autoregressive process with exogenous input, with {Z t } as its threshold process and s as its seasonal period. The vector of nonstructural parameters of the TSARX model is defined as . . . , (h (l) ) 2 ) and θ s = (k 1 , . . . , k l , K 1 , . . . , K l , q 1 , . . . , q l , r ′ , d, s, l) ′ is the vector of structural parameters of the model. θ = (θ x , θ z ) is denoted as the complete parameter vector, where θ x = (θ ′ ns , θ ′ s ) is the parameter vector in the TSARX model and θ z is the parameter vector of the Markov chain.
(ii) The conditional cumulative distribution function of X t given x 1,t−1 , z 1,t−1 , t > K, is given by where p t,j = P r(z t−d ∈ R j ), l j=1 p t,j = 1 and F t−1,j (· |· ) (j = 1, . . . , l) is given in (i). Now, being that {Z t } converges weakly to F , p t,j → p j = F (r j ) − F (r j−1 ) when t → ∞ for all j = 1, 2, . . . , l. Of course, if {Z t } is a strictly stationary process, then the cumulative distribution function of Z t is F for all t and, thus, p t,j = p j for all t and for all j = 1, . . . , l.
The theoretical considerations (i) and (ii) above are taken into account for the calculation of the generalized residuals of the TSARX model.
The model (1) can be considered as a special case of a TARX model, as described below. A TARX( l; k 1 + sK 1 , . . . , k l + sK l ; q 1 , . . . , q l ) model is of the form if r j−1 < Z t−d ≤ r j for each t ∈ Z, some j = 1, . . . , l, k j < s for all j, and {ϵ t } ∼ GW N (0, 1).
The complete parameter vector of the TARX model is conveniently defined as sKj ), for each j = 1, . . . , l, the TARX model given in (2) becomes the TSARX model given in (1). Related to the restrictions given regarding the coefficients of A j2 of a TARX model, in a later subsection a statistical test is presented to determine multiplicative seasonality by regimes

The Conditional Likelihood Function
Let x 1,T = (x 1 , . . . , x T ) and z 1,T = (z 1 , . . . , z T ) be observed data vectors for the processes {X t } and {Z t }, respectively, where T is the length of the sample period. Following Nieto (2005), it is assumed that the probabilistic mechanism that generates z 1,T does not depend on θ x and the joint distribution of x 1,T conditional on z 1,T and θ, does not depend on θ z . With these conditions it can be seen that, initially, it is possible to estimate the parameters of {Z t } and then conditional on these parameters, we proceed to estimate the parameters of the TSARX model considering the likelihood function 1 Our interest is to calculate the conditional likelihood function of a TSARX(l; k 1 , . . . , k l ; K 1 , . . . , K l ; q 1 , . . . , q l ) s model given in (1); that is to say, f (x 1,T | z 1,T , θ x ). Additionally, it is assumed that the first K values of x 1,T ; that is, x 1,K = (x 1 , . . . , x K ) are fixed and known with K = max{k j +sK j , q j , j = 1, . . . , l}; consequently, where . . , t nj ,j } are the indexes in time, where r j−1 < z t−d ≤ r j and n j is the number of observations in the regime j, j = 1, . . . , l. Henceforth, the conditional likelihood function of a TSARX model is given by (3).

Estimation of Parameters of a TSARX Model
Stochastic search ideas are used to identify autoregressive orders in a TSARX model with l regimes; in particular, the Gibbs variable selection method given by Dellaportas, Forster & Ntzoufras (2002) and denoted by GVS. This method achieves the estimation of autoregressive orders in a single stage and can involve these estimates in the MCMC sampling scheme, along with the remaining parameter estimates, assuming that the number of regimes is known. The GVS method is a hybrid between the method for the selection of variables proposed by Kuo & Mallick (1998), denoted KM, and the stochastic search for the selection of variables (SSVS) given by George & McCulloch (1993). These and other Bayesian model selection techniques can be reviewed in O'Hara & Sillanpä (2009). So & Chen (2003) use the SSVS method to determine the autoregressive orders and the remaining parameters, when the number of regimes in the SETAR models is known. Calderón & Nieto (2017) use the KM and GVS methods in multivariate TARX models.
Model (1) is rewritten with the vectors S j , M j and N j , j = 1, . . . , l as for some j = 1, . . . , l and the other specifications given in subsection 2.1 are met, except that the vector of structural parameters of the model is now To use the GVS method, it is necessary to assume that: where P r(s ij = 1) = π ij , P r(m uj = 1) = uj and P r(n vj = 1) = ω vj , i = 0, 1, . . . , k j , u = 1, . . . , K j , v = 1, . . . , q j , j = 1, . . . , l. According to George & McCulloch (1993), the interpretation of this formulation is as follows. First, we i would probably be so small that it could be safely estimated by 0. Second, we set 1/η 2 ij φ 2 ij large so that if s ij = 1, then a non-0 estimate of a (j) i should probably be included in the final model. A similar deduction can happen for the other mixtures of normal distributions. Additionally, those mixtures of normal distributions can be carried to a multivariate structure in the following way: The conditional likelihood function of a TSARX model with l regimes considering the vectors S j , M j and N j , j = 1, . . . , l is given in (3), except now, with, k j , K j , and q j , j = 1, . . . , l, defined previously.
Once MCMC samples are obtained for S j , M j , N j , j = 1, . . . , l, the autoregressive orders of the model can be identified considering a particular combination of S j , M j and N j , j = 1, . . . , l which possesses the highest posterior probability.
Prior independence is assumed for A j , B j and C j so V Aj0 = I kj +1 , V Bj0 = I Kj and V Cj0 = I qj , j = 1, . . . , l (I is the identity matrix). Also, it is assumed that each autoregressive coefficient is equally likely to be included in the model; then, π ij = 1/2, uj = 1/2 and ω vj = 1/2, The prior distributions for parameter groups A j , B j , C j , (h (j) ) 2 , j = 1, . . . , l, r and d are: (iii) For C j | N j a multivariate normal distribution with mean 0 qj and covariance (D Nj D Nj ) −1 , j = 1, . . . , l.
(v) r = (r 1 , . . . , r l−1 ) ′ follows the constant prior distribution f (r) = 1 C I(A) where C = . . . A dr 1 . . . dr l−1 , and A is a region that satisfies: (1) a ≤ r 1 < . . . < r l−1 ≤ b, with a and b convenient quantiles of Z t−d , and (2) each regime contains at least one H% from the sample of Z t−d (this to ensure valid inference). The value of H is preset.
(vi) d follows a discrete uniform distribution, with probability mass function Combining the conditional likelihood function for a TSARX model, considering the vectors S j , M j and N j , j = 1, . . . , l, with each of the prior distributions given in (i) to (vi), the posterior distributions of the unknown parameters of the model are obtained.
. . , l, r and d are generated from convenient complete conditional distributions as follows.

Initialize:
. , x tn j ,j ) ′ and n j the number of observations in the regime j, j = 1, . . . , l.
Here θ x−Aj is the parameter vector θ x without the subvector A j .
where f (· |· ) is the conditional likelihood function of the TSARX model taking into account, that S j , M j and N j , j = 1, . . . , l are latent indicator variables and A is the region specified in the prior distribution for r.
7. Generate d from the multinomial distribution, with probability mass function . . , l individually based on the Bernoulli distribution, with probability mass function

Generate vector
All the posterior distributions of Algorithm 1, except that of r, are standard. To simulate values of r, a Gaussian random walk Metropolis algorithm with step sizes σ r1 , . . . , σ r l−1 is used (Metropolis, Rosenbluth, Rosenbluth & Teller 1953). We set the MCMC sample size N sufficiently large, discarding the burn-in iterates, and keep the last G = N − M iterates for posterior analysis.
The means and standard deviations of the posterior distributions of A j , B j , C j , (h (j) ) 2 , j = 1, . . . , l and r are calculated to obtain point estimates and a measure of uncertainty for such estimates, respectively. Also, with the samples generated, 100(1 − α)% credible intervals or 100(1 − α)% Highest-Density Regions (HDR), for α, 0 < α < 1, can be obtained. The estimation of the delay parameter d is taken as the posterior mode and for the estimates of the vectors of zeros and ones, S j , M j , N j , j = 1, . . . , l the highest occurrence probabilities are taken.
In order to examine each of the parameters, if the MCMC samples obtained converge to the posterior distribution, the Geweke's Z-score plot (Geweke 1992) of the geweke.plot function of coda in R software (Plummer, Best, Cowles & Vines 2006) is used. Also, we review the evolution of the point estimate (posterior averages) of each of the parameters, when the number of iterations is increased, by means of the cumuplot function of coda in R. This plot presents the sample quantiles 0.025, 0.50 and 0.975, as a function of the number of iterations (accumulated averages) of the Markov chains.

A Bayesian Test of Multiplicative Seasonality by Regimes
De Gooijer & Vidiella-i Anguera (2003) using a frequentist approach, present a model that captures multiplicative seasonality by regimes and nonlinearity of thresholds simultaneously in SETAR models.
In this Subsection, to test the presence of multiplicative seasonality in each of the regimes of the TARX model, a Bayesian statistical test is proposed where an approximate technique is used for the calculation of the Bayes factor. In particular, we use the Savage-Dickey distributions ratio given in Dickey (1971) paper, which is an approximation of the Bayes factor (see, Verdinelli & Wasserman, 1995).
The methodology consists in comparing the TSARX model with a model where such seasonality is not captured (a TARX model). That is, that the parameter vector A j2 defined for the TARX model given in (2), equals the vector sKj ), for each j = 1, . . . , l, so that a TSARX model is a special case of a TARX model.
Formally, we want to test the hypotheses H 0j : A j2 = A j20 versus H aj : A j2 ̸ = A j20 , for each j = 1, . . . , l. The Bayes factor is given by where, respectively, f (·) and f (· | ·) denote the prior distribution and the conditional likelihood function of the TARX and θ x0 summarizes the parameters of the TARX model in the case A j2 = A j20 , for j = 1, . . . , l.
The Bayes factor (4) is equal to for j = 1, . . . , l; that is, the ratio of the weights of the conditional posterior distribution of A j2 and the weights of the prior distribution of A j2 , both evaluated in A j20 , j = 1, . . . , l.
The weight of the marginal posterior distribution of A j2 in A j20 can be obtained with the calculation of the average of the conditional posterior distribution of A j2 over the MCMC samples of the other parameters; thus, for j = 1, . . . , l, if the complete conditional posterior distribution of A j2 has a known standard distribution (see, Gelfand & Smith, 1990) and θ x−Aj2 is the parameter vector θ x without subvector A j2 . Here it is necessary to obtain MCMC samples of the posterior distributions of A j2 . Given the number of regimes l and autorregresive orders k j + sK j , q j , j = 1, . . . , l, the remaining parameters of the TARX model defined in (2), are estimated using a Gibbs sampler.
Returning to the Bayes factor calculation, because the parameters θ j2 and V j2 are used from the posterior distribution and θ j20 and V j20 from the prior distribution for A j2 , for the calculation of the Bayes factor by regimes as indicated in equation (4). Thus, , (strong or very strong evidence in favor of H 0j ) for each j = 1, . . . , l, according to the decision criterion given by Kass & Raftery (1995).

Estimation of the Number of Regimes
The identification of l, the number of regimes, can be seen as a problem of Bayesian model choice. The decision among models M j , j = 1, . . . , l 0 , M j being a TSARX model with j regimes and l 0 being a maximum number of regimes, is performed by the calculation of posterior probabilities, where the model M j is chosen if P r(M j | x K+1,T , z 1,T ) is the highest one over j = 1, . . . , l 0 . Congdon (2006) developed the following posterior probability estimator, which requires separate independent MCMC samples for each of the candidate models and is denoted as where φ is the prior distribution for model k and f (x K+1,T | x 1,K , z 1,T , θ (i) k ) is the conditional likelihood function for M k , for each k = 1, . . . , l 0 . For numerical efficiency, the calculation of (5) is based on scaled versions of the log-likelihoods. The scale uses the maximum L (i) max of the log-likelihoods in each iteration and subtracts this from each log-likelihood model, that is Some results support that the method is quite accurate and powerful for model selection. For example, Chen et al. (2010) compare TARX models with GARCH errors for different regimes, while Gerlach & Chen (2008) use it for STAR models with GARCH errors. The effectiveness of this approach is examined, comparing it with the deviance information criterion (DIC), which is described below. Spiegelhalter, Best, Carlin & Van Der Linde (2002) propose a Bayesian model comparison criterion defined as: where . . , G, j = 1, . . . , l 0 . The model M j , j = 1, . . . , l 0 is chosen if its DIC is the smallest among the candidate models.

Diagnostic Checking
Testing for model adequacy is an important part of any time series analysis. Based on MCMC techniques and Bayesian inference, it is natural to consider techniques that take into account the uncertainty of the parameter estimators in the calculation of the residuals. In this paper, the method of Gerlach, Carter & Kohn (1999) is used, which is based on MCMC techniques, importance sampling and the time series is the cumulative conditional distribution function given in Subsection 2.1, and K = max{k j + sK j , q j , j = 1, . . . , l}. Gerlach et al. (1999) show that for k ≥ t, converges in distribution to u t when G → ∞. Here θ (i) k is the i-th MCMC iterated sampled from the posterior distribution f (θ k | x 1,k , z 1,k ) for a given k and f (· | ·) is the conditional likelihood function for the TSARX model with l regimes. Because the variance of u t in equation (7) increases with (k − t), it is necessary to calculate u t with t reasonably close to k. This can be achieved by increasing k sequentially, say to 100, 200, . . . , T and evaluate u t using the equation (7) with (k − t) not greater than the increments. Based on the convergence properties of u t , the generalized residuals v t = Φ −1 ( u t ) (Φ is the standard cumulative normal distribution function) are approximately independent and identically distributed N (0, 1), under the correct model. Standard diagnostic tests can be applied to v t such as the Ljung-Box test to determine autocorrelation, the Jarque-Bera test to verify normality, the CUSUM plot to determine the correct specification of model and the CUSUMSQ plot to check marginal heteroscedasticity, among others.  and Chen et al. (2010) conduct the diagnostic check on TARX models with GARCH errors using generalized residuals and their results were very satisfactory for determining the adequacy of the model.

Forecasting
Here, the idea is to find X T +h = E(X T +h | x 1,T , z 1,T ), where h ≥ 1 is the forecast horizon, x 1,T = (x 1 , . . . , x T ), z 1,T = (z 1 , . . . , z T ), and T is the length of the sample period. X T +h is the best prediction in the sense of the minimum value that minimizes the quadratic lose function and its analytical expression is difficult or impossible to obtain in the context of nonlinear time series (see, . For this reason, we focus on the so-called predictive distribution of X T +h , that is, f (x T +h | x 1,T , z 1,T ). The predictive distributions of interest can be obtained marginally from the joint predictive distribution with θ x the parameters vector of the TSARX model and θ z the parameters vector of the variable Z t . Previous work in this direction is presented in Nieto (2008) and Vargas (2012), among others. It will be assumed in what follows, that the specifications given in previous subsections, with respect to the chosen model, are fulfilled. To obtain samples from the predictive distribution (8), we take into account that the joint predictive distribution can be expressed as is the posterior distribution of the parameters of the TSARX model; Also, x 1,T +i−1 = (x 1 , . . . , x T +i−1 ) and z 1,T +i−1 = (z 1 , . . . , z T +i−1 ). Thus, from the joint predictive distribution, z T +1 , x T +1 , z T +2 , x T +2 , . . . , z T +h , x T +h can be extracted sequentially as follows (Congdon 2007): Algorithm 2.

Generate a MCMC sample of θ (k)
x of the posterior distribution f (θ x | x 1,T , z 1,T ).

Generate z
x ).

Generate x
x ).
With the samples {z T +h }, k = 1, . . . , G, with G being the length of the chain and h ≥ 1, we calculate the means (point forecasts), the standard deviations, and credible intervals or HDR, for both variables.
The procedure described above includes the uncertainty of the estimation of the model parameters in the calculation of the forecasts, as is discussed by Vargas (2012) in TAR models.
In summary, our proposed methodology for fitting a TSARX model to a time series of interest that exhibits seasonality, consists of the following steps: Step 1. An exploratory analysis is made for the threshold and interest variables, to detect the presence of seasonality. For this, plots of the time series, correlograms, and plots of the time series for seasonal periods are used.
Step 2. Maximum autoregressive orders are fixed, which are obtained by fitting candidate ARX models to the data and choosing that order with minimum DIC. The threshold nonlinearity test proposed by Tsay (1998) is carried out and the presence of multiplicative seasonality is tested, using the statistical test proposed in Subsection 2.4.
Step 3. We proceed to estimate the number of regimes, autoregressive orders and the remaining parameters of the TSARX model, following the procedures described in Subsections 2.3 and 2.5.
Step 4. The estimated model is validated, using the generalized residuals as described in Subsection 2.6.
Step 5. Finally, the fitted model is used to obtain forecasts of the threshold and interest variables, as indicated in Subsection 2.7.
All procedures were implemented in R package (Team 2016).

A Simulated Example
In this section we present a simulation example in order to ilustrate the performance of the proposed methodology 2 .

The model
We assume that {X t } follows the TSARX(2; 1, 1; 2, 1; 1, 3) 12 model given by where {Z t } obeys the AR(1) model Z t = 1.80 + 0.60Z t−1 + a t , with {a t } a zeromean Gaussian white noise process with variance 1, which is independent of {ϵ t } and r = 4.46 is the median of Z t−2 . Notice that {Z t } is a homogeneous first-order Markov chain, with kernel density function f 1 (z t | z t−1 , θ z ) and corresponds to a N (1.80 + 0.60z t−1 , 1).

Exploratory Analysis of the Data
Time series plots, the sample autocorrelation function (ACF), the sample partial autocorrelation function (PACF), and the average mutual information index (AMI) were used to preliminarily determine the presence of seasonality.
The AMI measures the mutual dependence between variables W i and W i−τ , τ = 0, 1, . . . , τ * , where τ * is a specified maximum number of lags. If W i and W i−τ are independent, then knowing W i does not give information about W i−τ and viceversa, so their mutual information is zero. The mutual function of the package tseriesChaos in R (Di Narzo & Di Narzo 2013) is used to estimate the AMI.
In Figure 1, the two simulated time series are ploted, with a sample size T = 600 observations. Seasonal behavior is not detected at first sight. For {X t } (see Figure 2) it is observed that significant autocorrelations can be identified in the first lags and lag 12, essentially in the PACF. The AMI only shows the first three lags to be significant. For {Z t } (see Figure 2), only the first lag is significant, as indicated by the ACF, PACF and AMI plots. Now, looking at Figure 3

Estimation of the TSARX Model Parameters
100 simulated time series of sample size T = 600 are generated from model (9). For each simulated time series, the MCMC sampler is run for 12.000 iterations and the first 6.000 iterations are taken as the burn-in period. The posterior inference is based on the remaining 6.000 iterations. The prespecified maximum autoregressive orders k 1 , K 1 , and q 1 , are chosen from an ARX(k 1 + sK 1 ; q 1 ) model that best fits the simulated data, according to the minimum DIC. We take candidate ARX models, with k 1 = 0, 1, 2, 3, 4, 5, K 1 = 0, 1, 2, 3, and q 1 = 0, 1, 2, 3, . . . , 15 and s = 12. The values k 1 = 2, K 1 = 2, and q 1 = 3 were obtained. The Tsay threshold nonlinearity test (Tsay 1998) is applied for delay d = 0, 1, 2, 3. The null hypothesis of a linear ARX(24; 3) model is rejected for d = 2 at the 1% significance level. In the Bayesian test of multiplicative seasonality by regimes, it is required to estimate parameters of a TARX model with two regimes, with k 1 = k 2 = 2, K 1 = K 2 = 2, q 1 = q 2 = 3, and s = 12.  For the hyperparameters the following values were taken: θ j10 = 0 kj +Kj +1 ,  [(h (j) ) 2 ] (0) = 0.20, j = 1, 2; d (0) = 3; the step size to execute the Metropolis algorithm for drawing samples of r was σ r = 0.025, and r (0) the 0.50 quantile of Z t−d . When we run the Bayesian test for the 100 simulated series, the presence of multiplicative seasonality it is detected in both regimes, with average values of 2ln(F B 1 )=3725.30 and 2 ln(F B 2 ) = 15072.52. In the identification of the number of regimes, separate and independent MCMC samples are used for the TSARX candidate models, with l = 1, 2, 3, 4, 5 regimes (l 0 = 5). Additionally, we set φ ij = 25, η ij = 1.50 for all i, j, χ uj = 25, ξ uj = 1.50 for all u, j, ψ vj = 25, λ vj = 1.50 for all v, j, i = 0, 1, 2, u = 1, 2, v = 1, 2, 3, and j = 1, 2, 3, 4, 5. These values are the same as those chosen by Calderón & Nieto (2017). Furthermore, ν j and λ j , j = 1, 2, 3, 4, 5; a, b, H and d 0 are the same hyperparameters that were assumed for the TARX model used in the seasonality test. Now, the initial values for the parameters under the null hypothesis were A (0) [(h (j) ) 2 ] (0) = 0.20, j = 1, 2, 3, 4, 5; also, d (0) = 3, r (0) = q 0.50 {Z t−d } for two regimes, r (0) = (r 1 , With the step sizes σ r1 = σ r2 = σ r3 = σ r4 = 0.025 of the Metropolis algorithm we draw from the distribution of r. The results of the estimation of the number of regimes are shown in Table  1. Here we see the averages of the posterior probabilities given in (5), which are computed using each one of the 100 time series, and the average of the standard DICs given in (6), as an alternative to doing this task. We find that the true model is selected.
Considering l = 2 and maximum autoregressive orders k 1 = k 2 = 2, K 1 = K 2 = 2 and q 1 = q 2 = 3, the true model can be expressed as: = 1 1 0 1 1 1 0 0 1 1 0 1 0 1 1 1 , Table 2: The best three models identified using the GVS method. where the ones indicate that the coefficients associated with these latent variables are present in the model. The three best selected models, are reported in Table 2. Using the GVS method, the best model that is selected matches the true model with a posterior probability of 0.1801 and the second and third best models differ from the true model only in the presence or not of some coefficients of the exogenous input in the first regime. The three selected models are associated with d = 2. Also, the proportion of times that the model is declared as either the best or the second best are 88% and 95%, respectively.  (1) We proceed now, to estimate the other parameters of the TSARX model. The estimation of the autoregressive parameters in each regime, the variance weights by regimes, the threshold and delay values of the 100 replications are summarized in Table 3, showing each of the true values, plus the average, the standard deviation, the 0.025 and 0.975 quantiles of the 100 posterior means for each of the parameters. In addition, the percentage of coverage is given, which is the number of times that the true value of the parameter fall within the credible interval. For d, the average of the 100 posterior mode is given. All the averages of the posterior means are close to the true values and the 95% confidence intervals of the posterior means contain the true values. The percentages of coverage are reasonable and very close to 95%. Also the delay d = 2 is correctly identified with the sampling scheme considered.
In order to ensure that for each parameter, the MCMC samples converge to the posterior distribution, cumuplots function and Geweke's Z-score plots are used in one of the one hundred simulated time series. Because there are several parameters in the model, only some of them are shown (Figure 4), where we observe that convergence is reached in each case.

Model Diagnostics
For each one of the 100 simulated time series we estimate u t , t = 101, . . . , 600 (see equation 7). This is done by setting the value of k sequentially at 150, 200, 250, . . . , 600, and for each t such that k − t < 50; then, v t is calculated, which is the generalized residual at time t. The Ljung-Box tests for lags 1, 12 and 24 and the Jarque-Bera test are used to check whether v t violates the hypothesis that it is independent and identically distributed N (0, 1). We reports the number of nonrejections out of 100 replications (not shown in table). The model (9) is correctly fitted, it is good to see that the empirical sizes of the Ljung-Box tests are 96%, 96% and 92% and the Jarque-Bera test is 94%, respectevely, match the nominal values of 95% very closely. Using one of the 100 simulated time series, CUSUM and CUSUMSQ plots of the generalized residuals are presented in Figure 5, indicating that there is no evidence of an incorrect specification of the model nor the presence of marginal heteroscedasticity in the generalized residuals, respectively.

Forecasts Computation
Now, the forecast procedure developed in Subsection 2.7 is illustrated. In order to check the forecasting performance of the model, we considered as sample period 1-587 and the forecast horizon was 588-600. In Table 4, the simulated values, point forecasts, standard deviations and 95% credible interval averages for the variables X and Z are presented. It is worth noticing that these estimates are obtained via the averages of the corresponding statistics in each one od the 100 simulated time series. We put h = 1, 2, 3, 11, 12, 13. Additionally, in the last column of the table, the coverage percentage is presented.
We find that the forecasts are very close to the simulated values, all the simulated values fall within the corresponding 95% credible intervals, the standard deviations increase as the forecast horizon increases, and the coverage percentage is close to the 95% confidence level. The posterior distributions of the predictive distributions present a high variability. In this way, We assume that it is due to the use of priori noninformative distributions. Similar results to those found in Table  4, are given in Nieto (2008) and Vargas (2012), for TAR models. The behavior of the uncertainty measure of the long-term forecasts for TSARX models, maybe investigated, since it escapes the purpose of this paper.

An Empirical Application
In this section, we apply the proposed methodology to Colombian monthly economic time series. The variables are the total monthly unemployment rate, TD t , and a monthly coincident index that traces the Colombian economy, ISE t . In our analysis, we focus on the growth of the total monthly unemployment rate, X t = T D t − T D t−1 , as the interest variable, and the monthly growth of the logarithmic transformation of ISE t , Z t = (ln(ISE t ) − ln(ISE t−1 )) × 100 (in percent) as the threshold variable. The sample period is from February, 2000to December, 2016. The data is provided by the DANE, the Colombian official agency for the study and diffusion of official statistics.

Exploratory Analysis of the Data
The time series {X t } and {Z t } are given in Figure 6(a) and (b), respectively. {X t } presents annual positive peaks (Januaries), while {Zt} presents very pronounced annual negative peaks (Januaries). This indicates some seasonal behavior in both series, inherited from the original time series. It can also be observed that the series have no trends, but present seasonal cycles. In Figure 7, autocorrelograms and the AMI for X t and Z t are shown. In both cases, the samples ACF show seasonal lags that do not decrease and the PACF and AMI identify significant correlations in the first lags and lags 12, 24 and 36. Figure 8 suggestes fitting a multiplicative AR(26) model to the two time series, following the minimum BIC. Looking at the box plots, X t presents a distribution with higher values in the months of January, and a distribution with lower values in the months of March. On the other hand, for Z t there is a distribution with lower values in the months of January, and a distribution with higher values in the months of October. In this way, there is some empirical evidence about the presence of seasonality in the time series.  An adequate linear AR model fitted to the threshold process is

Estimation of the TSARX Model Parameters
To obtain the maximum autoregressive orders, the same candidate ARX models of Subsection 3.3 are considered. The values k 1 = 3, K 1 = 2 and q 1 = 2 were achieved. These values are fixed to apply the Tsay threshold nonlinearity test (Tsay 1998), the proposed Bayesian seasonality statistical test, the Congdon method and the GVS technique.
The MCMC sampler for posterior analysis is again based on 6.000 iterations after a burn-in period of 6.000 iterations. The hyperparameters and initial values for each of the parameters of the TARX models, are the same as those in Subsection 3.3. In Table 6, the results of the seasonality test are listed. A multiplicative TSARX model with two regimes is selected. For the case of l = 3 regimes, the second regime does not show the presence of multiplicative seasonality.  Now, we estimate the number of regimes and the autoregressive orders of a TSARX model. The hyperparameters and initial values of the parameters for each possible TSARX model are the same, as those described in Subsection 3.3. The posterior probability of the number of regimes of a TSARX with j = 1, 2, 3, 4, 5 regimes, is given in Table 7, where it can be seen that the number of regimes associated with the highest posterior probability is two. This result is confirmed with the smallest DIC value. In Table 8, the three best models are observed based on the highest posterior probabilities, for the latent variable matrices S j , M j and N j , j = 1, 2. The probabilities of the best three TSARX models with two regimes are 0.2122, 0.1536 and 0.1078. The three models are associated with d = 3 and they look very similar. Indeed, they have the same structure in the second regime and differ mainly in the entries of coefficients a (1) 3 and b (1) 2 . The best model is chosen and the remaining parameters are estimated, which are shown in Table 9. Following the obtained results, the TSARX model is explicitly given by where the model for Z t is given in model (10).  An interpretation of the fitted TSARX model is the following: first of all, only two regimes for the growth of the Colombian total monthly unemployment rate are detected, which could be termed as regimes of the monthly growth of the logarithmic transformation of the ISE t , low and high. In the low regime, the most pronounced negative peaks are also included. In addition, it is interesting to note that for the first regime, the intercept is negative and affects the epochs of recession of the economy by 22%, this reflects the economic impact that the growth of the unemployment rate has in that period; on the contrary, in the periods of expansion of the economy that are modeled by the second regime, the intercept is zero, giving evidence of a rapid economic recovery. Another interesting characteristics are the following: (1) the multiplicative AR model of the second regime has autoregressive orders higher than those of the first regime, (2) the coefficients of the nonseasonal component have negative signs, and the coefficients of the seasonal component have positive signs of similar magnitude, in the two regimes of the model, (3) the variances of the regimes are similar, and (4) the monthly growth of the logarithmic transformation of the index, affects the growth of the unemployment rate, which took place three months earlier. Table 9: Estimation of parameters for empirical data.
In Figure 9 the cumuplots function and Geweke's Z-score plots are shown, respectively, for some of the MCMC chains of the estimated parameters. The other chains can be provided by the authors upon request. The convergence, to the posterior distributions, is almost immediate and clearly reached, after the burn-in period in each case considered.

Model Diagnostics
In Table 10 we present the results of applying the Ljung-Box test for lags 1, 2, 3, 12 and 24 and the Jarque-Bera statistical test to the generalized residuals. In all cases, the null hypotheses of zero autocorrelations and normality are not rejected, under the 5% significance level. These results suggest that the TSARX model with two regimes fit reasonably well the empirical data. Figure 10 shows the CUSUM and CUSUMSQ plots of the generalized residuals, which complements the good fit of the model.

Forecasts of the Variables
As in the simulation example, we compute now forecasts for the variables X and Z. Our goal is double: first, to compute point predictions and second, to check the model performance for forecasting. We use 6.000 MCMC iterations in obtaining samples from the corresponding predictive distributions. We present the results in Table 11, beginning in January, 2017 (h = 1) until January, 2018 (h = 13), period for which we include additional observations. The table includes the observed values, point forecasts, standard deviations of the predictive distributions and 95% confidence-level credible intervals. It can be noted that in both variables, (1) the point forecasts are very close to the observed values, (2) the standard deviations increase with larger values of h, and (3) all the observed values fall within the credible intervals. In future research, we will compare forecasts from the TSARX models with those from linear seasonal models and other nonlinear seasonal models, following the ideas given in Franses & Van Dijk (2005) and Vaca (2018).

Conclusions
In this paper, an open-loop threshold autoregressive model is specified for explaining the presence of seasonality in a nonlinear time series. The autoregressive model is, essentially, multiplicative but depending on the regimes. A Bayesian test is proposed for detecting this kind of seasonality. To fit the model to a time series, a full Bayesian methodology is developed that consists in the phases of (i) identifying the model, (ii) estimating the unknown parameters, (iii) checking the diagnostic via generalized residuals, and (iv) using the model to compute forecasts of the variables. The proposed methodology was illustrated with a simulated example and a real-data application. In future research, a more extensive characterization of the model may be conducted, in order to detect in the data more closely this kind of seasonality. Also, an extension to nonlinear multivariate time series may be investigated.