Some Dynamic and Steady-State Properties of Threshold Auto-Regressions with Applications to Stationarity and Local Explosivity

Abstract

The purpose of this paper is to investigate the dynamics and steady-state properties of threshold autoregressive models with exogenous states that follow Markovian processes. Markovian processes are widely used in applied economics although their statistical properties have not been explored in detail. We use characteristic functions to carry out the analysis, and this allows us to describe limiting distributions for processes not considered in the literature previously. We also calculate analytical expressions for some moments. Furthermore, we see that we can have locally explosive processes that are explosive in one regime whilst being strongly stationary overall. This is explored through simulation analysis, where we also show how the distribution changes when the explosive state becomes more frequent although the overall process remains stationary. In doing so, we are able to relate our analysis to asset prices which exhibit similar distributional properties.

1. Introduction

The purpose of this paper is to investigate the dynamics and steady-state properties of threshold autoregressive models with exogenous states that follow Markovian processes. These models fall within the class of regime-switching models, which have become increasingly popular in applied economics and finance. Initially introduced by Goldfeld and Quandt (1973) and Tong (1978), regime-switching models have been used in economics and finance for a wide variety of applications including forecasting exchange rates (Engel 1994), understanding price transmission (Goodwin and Harper 2000), detecting bubbles in the art market (Knight et al. 2014), and providing a metric of market efficiency (Ahmed and Satchell 2018). Hansen (2011) provides a concise summary of threshold autoregressive processes and their applications.

Hamilton (1989, 1990, 2010) has made seminal contributions to the theory and application of regime-switching models. As outlined above, this article discusses a particular class of regime-switching models. The problems we discuss appear to have much in common with Markov switching models, and Timmermann (2000) has provided a detailed analysis of moments and autocorrelations, which would include our model as “MSIII” in his terminology. However, his analysis does not address non-moment distributional properties or the non-existence of moment-generating functions (mgfs).

Indeed, Timmermann states in Appendix 1, page 103, that “The expressions for the cases where ${\epsilon }_t$ follows a t-distribution or a normal distribution are based on the moment-generating distributions for these distributions”, and this is confusing as it is known that the t-distribution does not have a moment-generating function. We therefore re-examine this model, allowing for the non-existence of moment-generating functions, and use the characteristic function (which will always exist) to derive various properties of the model.

Whilst we could carry out a similar analysis for the other models described in Timmermann (2000), our focus is on threshold auto-regression and the elusive search for explicit steady-state distributions. Prior to this article, Gonzalo and Gonzalez-Rozada (1997) described statistical properties of Threshold Autoregressive models of order 1 (TAR(1)) models. However, their analysis was restricted to a mixture of stationary and unit roots, whereas our analysis considers non-stationary roots as well. Another important contribution relevant to our work is Pourahmadi (1988). In Theorem 3.1 in his article, Pourahmadi (ibid.) discussed the covariance stationarity of a process similar to the threshold autoregressive process we consider in Section 4, in that he analyzed processes with a unit root and a zero root. However, the results derived in Section 2 below and the processes considered in Section 4 are applicable to a more general setting and we do not restrict ourselves to covariance stationarity. We acknowledge Pourahmadi’s contributions while deriving our results in Section 2.

Caner and Hansen (2001) and Kapetanios and Shin (2006) also considered similar processes, but their objective was to derive the distribution of unit root test statistics in the threshold autoregressive framework rather than the distribution of the underlying process. Our results build upon the results of Knight and Satchell (2011) and Ahmed and Satchell (2018): Both articles discussed theoretical moments for threshold autoregressive models with exogenous triggers. While Ahmed and Satchell (2018) only considered moments when the exogenous variable is independently and identically distributed, Knight and Satchell (2011) also considered a Markovian exogenous variable.

In addition to deriving theoretical moments, we also use simulation analysis to show how the distributions of threshold autoregressive models change when the process’s two states consist of one stationary state and an explosive state. Our analysis focuses on these models for two reasons. Firstly, we are able to derive a characteristic function for this case, thereby adding to the literature on analytical results for threshold autoregressive models. We believe this is a significant contribution to the statistical literature. Secondly, this class of models is of interest in the financial literature concerned with explosive roots. In particular, Theorem 6 in Section 2 specifically refers to a special case that is of interest to researchers working on asset pricing with some non-stationarity. The latter is what we consider for our simulation analysis, and thus we consider it the most important contribution of the current article. We also believe that these models can prove to be useful in the applied macroeconomic literature. The simulation analysis presented in this article will help the reader appreciate how these models can be useful in practice.

In applied Macroeconomics, for instance, Dynamic Stochastic General Equilibrium (DSGE) models often model shocks as AR(1) processes (see Schmitt-Grohe and Uribe 2004). The literature came under particular scrutiny after the financial crisis for its inability to simulate and thereby predict conditions and outcomes that were observed during the crisis. In addition to the absence of financial markets, such models are also restricted by their reliance on a stationary AR(1) model as a shock process, as these processes can rarely be used to study the kind of macroeconomic shocks that led to the financial crisis. On the other hand, these models will not have analytical or numerical solutions if the process is non-stationary.

We postulate that using a TAR(1) shock process which is stationary but nevertheless can exhibit locally non-stationary behavior can improve these models. Our work will enable calculation of moments for such shocks (where such moments exist), allowing the user to work with analytical solutions, or, if the user is deriving a numerical solution, ensuring that such a solution will exist. Indeed, some work has already started relying on Markov-switching DSGE models (Foerster et al. 2016). This paper complements the proposed methodology by enabling researchers to control and simulate shocks of specific variances.

There are further applications in the finance literature where TAR models have become popular. The applications may extend to forecasting oil prices through threshold models or in modelling exchange rate fluctuations. There are many areas where a TAR model and the characteristic functions we derive can provide more depth to the underlying analysis. To use one recent example, Aleem and Lahiani (2014) estimated a TAR model of exchange-rate pass through for Mexico. Their analysis was limited to an estimation of the threshold above which the pass through is greater. The characteristic functions from our article would have allowed them to estimate the volatility of exchange rates in their model, improving both their model and the resulting predictions. Similarly, in Ahmed and Satchell (2018), the empirical application relied on a Markovian exogenous trigger. The results from this article would have allowed them to derive moments of their empirical TAR(1). Corollary 2 in Section 2 offers one example of how the results of this article may contribute to applied and empirical work in finance.

The rest of the article is organized as follows. In Section 2, we present the derivation and formulae for characteristic functions of threshold autoregressive (1) models with exogenous Markovian triggers. Section 3 outlines the simulation methodology. A separate section is necessary since obtaining a sample from the steady-state distribution of a TAR model with a Markov-switching exogenous trigger is a non-trivial exercise. Section 4 presents and discusses simulation results and Section 5 concludes.

2. Moment-Generating Functions of TAR(1) Models with Exogenous Markov-Triggers

In this section, we introduce the threshold autoregressive model with a Markov-switching exogenous trigger. After introducing the model, we derive the moment-generating functions for this model and present some interesting results. We shift to characteristic functions after Theorem 3, where we do not have to assume the existence of all moments so that moment-generating functions for such processes need not exist. Characteristic functions on the other hand will always exist (Stuart and Ord 1994, chps. 3 and 4). Such processes are often used to model prices (particularly in finance), therefore, we refer to our model as a price process indicated by $p_t$.

The price process has a switching AR(1) form:

$$p_t={\alpha }_{t-1}+{\beta }_{t-1}{p}_{t-1}+{\sigma }_{t-1}{\eta }_t$$

(1)

where ${\alpha }_{t-1}$ is a switching drift, ${\beta }_{t-1}$ is a switching coefficient term, and ${\sigma }_{t-1}$ is a switching variance term for the error process. Let:

$${\alpha }_{t-1}=\langle \alpha ,z_{t-1}\rangle {\beta }_{t-1}=\langle \beta ,z_{t-1}\rangle {\sigma }_{t-1}=\langle \sigma ,z_{t-1}\rangle$$

where $\alpha$ is a vector containing all drift terms, $\beta$ is a vector of coefficients, and $\sigma$ is a vector of error standard deviations. In general, all the above vectors are k × 1, but we illustrate them when k = 3 for notational convenience. We assume that $z_t$ is Markovian and follows a multinomial distribution, and that ${\eta }_t$ has a moment-generating function $\psi \left(u\right)$ which is assumed to be location-scale. $z_t$ determines what state ${\alpha }_{t-1}$ and ${\beta }_{t-1}$ are in. In particular,

$$\begin{gathered} \alpha ={\left({\alpha }_1,{\alpha }_2,{\alpha }_3\right)}^{\prime } \beta ={\left({\beta }_1,{\beta }_2,{\beta }_3\right)}^{\prime } \sigma ={\left({\sigma }_1,{\sigma }_2,{\sigma }_3\right)}^{\prime } {\eta }_t~iid\left(0,1\right) \\ {z}_t\in \left\{e_1,e_2,e_3\right\} e_1=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right) e_2=\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right) e_3=\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right) \end{gathered}$$

$E\left[z_t|z_{t-1}\right]=P{z}_{t-1}$, where P represents the Transition Matrix for the Markovian state variable $z_t$. For econometric purposes, we envisage an exogenous continuous random variable $Z_t$ and constants ${\vartheta }_{0,}\dots .{\vartheta }_3$, so that $z_t=e_{j } \mathrm{if} {\vartheta }_{j-1,}\le Z_t <{\vartheta }_j$, i.e., when the continuous random variable $Z_t$ is between thresholds ${\vartheta }_{j-1}$ and ${\vartheta }_j$, and the value of the Markovian variable $z_t$ is equal to $e_j$. The nature of $Z_t$ determines $P$ and the kind of regime-switching model $Z_t$ is.

Here, $P=\left(\begin{array}{ccc}{p}_{11}& p_{12}& p_{13}\\ p_{21}& p_{22}& p_{23}\\ p_{31}& p_{32}& p_{33}\end{array}\right)$ is the transition matrix, which describes the probability of switching states. Note that $i^{\prime }P=i\prime$, where i is a vector of ones, and $P\pi =\pi$ where $\pi$ is the vector of stationary (steady-state) probabilities and $i^{\prime }\pi =1$.

$$p_{ji}=P\left(z_{t+1}=e_j|z_t=e_i\right)=P\left(z_1=e_j|z_0=e_i\right) 1\le i,j\le 3$$

so that the Markov Chain is stationary. Whilst we can estimate P by counting frequencies, we can also hypothesize a Markov process for $Z_t$ and then integrate over the appropriate rectangle of the probability density function of ($Z_t,Z_{t+1})$.

We now consider $\exp \left(u{p}_t\right)$ in order to derive the moment-generating function for $p_t$. Here, $u\in ℝ$.

The moment-generating function of $p_t$ is defined by ${\varphi }_t\left(u\right)=$ $E\left[\exp \left(u{p}_t\right)\right]$. Our aim is to determine a recursion for ${\varphi }_t\left(u\right)\in R$.

Now,

$$z_t=P{z}_{t-1}+{\nu }_t\in R^3$$

(2)

where,

$$E\left[{\nu }_t|z_{t-1}\right]=\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right)=0\in R^3.$$

From (1), we have that:

$$exp \left(u{p}_t\right)=exp\left[u\langle \alpha , z_{t-1}\rangle \right] exp\left[u \langle \beta , z_{t-1}\rangle p_{t-1}\right]exp\left(u\langle \sigma , z_{t-1}\rangle {\eta }_t\right)$$

Using iterated expectations, we find that the moment-generating function for $p_t$ is:

$${\varphi }_t\left(u\right)=E(E(exp\left[u\langle \alpha , z_{t-1}\rangle \right] exp\left[u \langle \beta , z_{t-1}\rangle p_{t-1}\right]exp(u\langle \sigma , z_{t-1}\rangle {\eta }_t|z_{t-1}))$$

For functions F ($p_{t-1}, Z_{t-1})$, we note that:

$$E(F(p_{t-1}, Z_{t-1}\left)\right) = E({{\displaystyle \sum }}_{j=1}^{k}F(p_{t-1}, e_j) {\pi }_j)$$

by the law of total probability.

Thus,

$${\varphi }_t\left(u\right)={{\displaystyle \sum }}_{j=1}^k{\pi }_j(exp\left[u\langle \alpha , e_j\rangle \right] {\varphi }_{t-1}\left[\langle \beta , e_j\rangle u\right]\psi \left(\langle \sigma , e_j\rangle u) \right)$$

(3)

is a dynamic recursion for the moment-generating function of $p_t.$

Steady-State Distribution under Markovian States

The above discussion leads to the following result.

Theorem 1.

Assuming a steady state for prices, denote $E\left(exp\left(u{p}_t\right) \right) = \varphi \left(u\right)$ and $E\left(exp\left( u{\eta }_t\right) \right) = \psi \left(u\right)$ as the appropriate mgfs (or characteristic functions with a trivial definitional change). Then,

$$\varphi \left(u\right)={{\displaystyle \sum }}_{j=1}^{k}exp \left(u{\alpha }_j\right) \varphi \left({\beta }_{j}u\right) {\pi }_j \psi \left({\sigma }_{j}u\right)$$

(4)

is the steady-state relationship.

We can use Theorem 1 to arrive at analytical expressions for different moments of the process, $p_t$.

Define:

$$\begin{gathered} {\mu }_B={{\displaystyle \sum }}_{j=1}^k {\pi }_j{\beta }_j; {\mu }_{2B}={{\displaystyle \sum }}_{j=1}^k {\pi }_j{\beta }_j^2; {\sigma }_B^2={\mu }_{2B}-{({\mu }_B)}^2 \\ {\mu }_{\alpha }={{\displaystyle \sum }}_{j=1}^k {\pi }_j{\alpha }_j; {\mu }_{\sigma }={{\displaystyle \sum }}_{j=1}^k {\pi }_j{\sigma }_j; {\mu }_{\alpha B}={{\displaystyle \sum }}_{j=1}^k {\pi }_j{\beta }_j{\alpha }_j\mathrm{etc}. \end{gathered}$$

We differentiate (4) once to obtain the first moment of $p_t$, and we get:

$$E\left(p_t\right)=\frac{{\mu }_{\alpha }+{\mu }_{\sigma }E\left({\eta }_t\right)}{1-{\mu }_B}$$

Differentiating the mgf a second time gives us the second moment of $p_t$:

$$E\left(p_t^2\right)=\frac{{\mu }_{2\alpha }+2{\mu }_{\alpha B}E\left(p_t\right)+2{\mu }_{\sigma B}E\left({\eta }_t\right)E\left(p_t\right)+2{\mu }_{\alpha \sigma }E\left({\eta }_t\right)+{\mu }_{2\sigma }E\left({\eta }_t^2\right)}{1-{\mu }_{2B}}.$$

Further calculations and simplifications lead to an expression for the variance of $p_t$:

$$\begin{gathered} Var\left(p_t\right) = E(p_t^2)-{\left(E\left(p_t\right)\right)}^2 \\ Var\left(p_t\right)=\frac{{\mu }_{2\alpha }+2{\mu }_{\alpha B}E\left(p_t\right)+2{\mu }_{\sigma B}E\left({\eta }_t\right)E\left(p_t\right)+2{\mu }_{\alpha \sigma }E\left({\eta }_t\right)+{\mu }_{2\sigma }E\left({\eta }_t^2\right)}{1-{\mu }_{2B}} -{\left(\frac{{\mu }_{\alpha }+{\mu }_{\sigma }E\left({\eta }_t\right)}{1-{\mu }_B}\right)}^2 \end{gathered}$$

Various interesting results can be derived from Theorem 1 for plausible parameter values. We list one case below, but other results can be regarded as special cases. Here, we are concerned with the case where k = 2, ${\alpha }_j=0, {\beta }_1=\beta$, ${\beta }_2=0, \mathrm{and}{\sigma }_2=0$, so that the mgf function becomes:

$$\varphi \left(u\right)=\varphi \left(\beta u\right)\left(\pi +\left(1-\pi \right)\psi \left(u\right)\right)$$

(5)

Corollary 1.

If $\psi \left(u\right)$ is the moment-generating function of a negative exponential with parameter $\lambda ,$ and ${\alpha }_j=0,{\sigma }_1=\sigma , {\sigma }_2=0$, and $\pi =\beta$ where ${\beta }_1={\beta }_2=\beta$ which is less than 1, then $\varphi \left(u\right)=\frac{\lambda }{\lambda -u}$, i.e., a negative exponential random variable with parameter $\lambda$.

Proof. 

To show that Equation (5) has a solution for some $\psi \left(u\right)$, we consider the negative exponential function, i.e., we assume that the disturbance term is distributed as a negative exponential with parameter $\lambda$. The corresponding moment-generating function, $\varphi \left(u\right),$ for this disturbance term is $\frac{\lambda }{\lambda -u}$ if we further assume that $\beta =\pi$ and that $0\le \beta <1$. Note that Equation (5) corresponds to the situation where the $\beta$ coefficient does not move across states, but the standard deviation of the disturbance term does, i.e., ${\sigma }_1=\sigma$ and ${\sigma }_2=0$. Our result implies that $\sigma =\frac{1}{\lambda }$. □

For our distributional assumption regarding the disturbance, the corresponding moment-generating function is (taking into consideration the two states):

$$\left(\pi +\left(1-\pi \right)\psi \left(u\right)\right)=\pi +\frac{\left(1-\pi \right)\lambda }{\lambda -u}$$

Substituting this in (5) and using the trial solution, $\varphi \left(u\right)=\frac{\lambda }{\lambda -u},$ we have:

$$\frac{\lambda }{\lambda -u}=\left(\frac{\lambda }{\lambda -\beta u}\right)\left(\pi +\frac{\left(1-\pi \right)\lambda }{\lambda -u}\right)$$

If we further assume that $\beta =\pi$, the LHS and RHS are equal, thereby proving our result.

We recognize the solution as being a negative exponential auto-regression of degree 1 (NEAR(1)). These models were investigated in detail by Gaver and Lewis (1980), which also included earlier references to related models. We note that the same arguments could be applied to Gamma random variables with integer degrees of freedom.

The attractiveness of these models is that they are AR(1) models where the underlying process is always positive, and hence, can be used to model equity or bond prices in finance. Our version is a slight extension of existing NEAR(1) models, in that Corollary 2 will be consistent with a Markov process for the state process rather than an i.i.d. process, as in the current NEAR(1) literature.

Furthermore, if $Z_{t-1}$ is Markovian with transition matrix P such that $\pi =P\pi ,$ then:

$${\pi }_j= {\displaystyle \sum }_{m=1}^k{\pi }_m{P}_{jm},$$

and

$$\varphi \left(u\right)= {{\displaystyle \sum }}_{j=1}^k\exp \left(u{\alpha }_j\right) \varphi ({\beta }_{j}u) \psi \left({\sigma }_{j}u\right) {{\displaystyle \sum }}_{m=1}^k{P}_{jm} {\pi }_m$$

Alternatively,

$$\varphi \left(u\right)=\langle \exp \left(u\alpha \right)\varphi \left(\beta u\right) \psi \left(\sigma u\right), P\pi \rangle$$

(6)

There are a number of observations relevant to (4) and (6) which we present below:

Theorem 2.

Since $\pi = P\pi$ has multiple solutions for $P$ given $\pi$, these different $P\prime s$ do not change the solution to Equation (6). As an example, for $k=2$, suppose $\pi =0.5$. It then follows that $P_{11}=P_{22}$, but if $P_{11}=0.2 or 0.8$ in this context, the steady-state distribution will be unaffected except through a change in position.

The steady-state values are equal (i.e., 0.5), and thus, such changes in the structure of the transition matrix should not influence steady-state values. Note, however, that this does not say anything about the speed at which the two processes in this example converge to the steady state. For more on speed of convergence, refer to Rosenthal (1995). Since the processes converge to the steady state through different paths, simulating the steady state becomes a non-trivial procedure, as explained in Section 3.

Theorem 3.

Suppose that in Equation (6), $\alpha$ is zero, and $\varphi \left(u\right) and \psi \left(u\right)$ are infinitely differentiable moment-generating functions and that the variance of the error process is constant.

Theorem 4.

$If \varphi \left(u\right) is symmetric$, then $\psi \left(u\right)$ is symmetric. The proof is trivial.

Theorem 5.

If $\psi \left(u\right)$ is symmetric, then $\varphi \left(u\right)$ is symmetric (proof by induction on Taylors series terms). We shall prove that all odd moments are zero:

Proof. 

$$\varphi \left(u\right)=\left({\displaystyle \sum }_{j=1}^k \varphi \left({\beta }_{j}u\right) {\pi }_j \psi \left(u\right)\right)$$

The coefficient of $u^n$ for $\varphi \left(u\right)$,

$${\varphi }_n={{\displaystyle \sum }}_{j=1}^k{\pi }_j{{\displaystyle \sum }}_{s=0}^n{\varphi }_{n-s}{\beta }_j^{n-s}{\psi }_s$$

(7)

${\psi }_1=0$ implies that ${\varphi }_1=0$. We now suppose that ${\psi }_{2j+1}=0$, which implies that ${\varphi }_{2j+1}=0$ for $j=0,\dots ,k$ and consider ${\varphi }_{2j+3}$. From (7), the inductive hypothesis and the properties of products of odd and even numbers, the result follows. □

Theorem 6.

Suppose that ${\alpha }_j=0$, the variance of the error term is constant, and that we treat

$$\varphi \left(u\right)=\left({\displaystyle \sum }_{j=1}^k\varphi \left({\beta }_{j}u\right){\pi }_j \psi \left(u\right)\right)$$

(8)

as a statement about characteristic functions. Then, if at least one of the${\beta }_j\prime s$is greater than 1 and all of them are non-negative, for some$n$, the$nth$moment will not exist. The proof follows from using (7) again and noting that${\varphi }_n$(the nth differential of$\varphi \left(u\right)$), which is proportional to the$nth$moment (if it exists), can be expressed as:

$${\varphi }_n\left(1-{{\displaystyle \sum }}_{j=1}^k{\pi }_j{\beta }_j^n\right)={{\displaystyle \sum }}_{j=1}^k{\pi }_j{{\displaystyle \sum }}_{s=1}^n{\varphi }_{n-s}{\beta }_j^{n-s}{\psi }_s$$

The requirement for the existence of ${\varphi }_n$ is that ${{\displaystyle \sum }}_{j=1}^k{\pi }_j{\beta }_j^n<1$, which cannot hold for a large enough n under the assumptions of Theorem 6. This result links local explosivity to fat tails. Thus, processes with locally explosive states will cease to have moments once $n$ becomes sufficiently large.

Ahmed and Satchell (2018) derived similar conditions for the existence of a mean and variance for a TAR(1) process with an independently distributed exogenous trigger for state switching. We have generalized the result for the $n\mathrm{th}$ moment and for an exogenous trigger that is Markovian. Pourahmadi (1988) arrived at a similar result in Section 2 of his article (see Equations (2.3)–(2.5) in Pourahmadi (ibid.)), but he carried out his derivation in the context of doubly stochastic processes, as opposed to the specific case of a threshold autoregressive process that we consider. Secondly, Pourahmadi was mainly interested in second-order stationarity, while we present results for the existence of all moments. Thus, we substantially improve upon the results contained in Pourahmadi (1988) and Ahmed and Satchell (2018). Below, we consider a special case which will be of particular interest to finance practitioners.

Corollary 2.

Assume that${\alpha }_j=0,and$the variance of the error term is constant, i.e.,${\sigma }_1={\sigma }_2=\sigma$.

Consider now the special case $k=2, {\beta }_1=1, \mathrm{and} {\beta }_2=0$. This is an important special case as it gives us a random walk in one regime and white noise in the other. Substituting into (8), we see that:

$$\varphi \left(u\right)=\left(\varphi \left(u\right)\pi +1-\pi \right)\psi \left(u\right)$$

This can be re-arranged to yield:

$$\varphi \left(u\right)=\frac{\left(1-\pi \right)\psi \left(u\right)}{\left(1-\pi \psi \left(u\right)\right)}$$

Since $\left|\psi \left(u\right)\right|\le 1$, $\pi \left|\psi \left(u\right)\right|$ < 1 and $\varphi \left(u\right)$ can be represented in terms of a valid series expansion which can be analyzed term by term. Indeed,

$$\varphi \left(u\right)=\left(1-\pi \right){\displaystyle \sum }_{j=0}^{\infty }{\pi }^j\psi {\left(u\right)}^{j+1}$$

(9)

The right-hand side is uniformly and absolutely convergent because of the Weirstrass M-test, and thus we can integrate term by term. Pourahmadi (1988) also considered this as a special case in his article and derives conditions for the existence of a variance and covariance (see Section 3 and Section 4 in his article).

We can now consider different choices for $\psi \left(u\right)$.

Suppose we have a normally distributed error term with a mean of zero and a variance of ${\sigma }^2$, i.e., $\psi \left(u\right)=\exp \left(-\frac{{\sigma }^2{u}^2}{2}\right)$. Then, $\psi {\left(u\right)}^{j+1}$ represents a normal random variable with a mean of zero and a variance of (j + 1) ${\sigma }^2$. We can identify the distribution of $p_t$ as an infinite weighted sum of normal random variables of increasing variances, but whose relative importance declines with a power of $\pi .$ This process was analyzed in Knight and Satchell (2011) and extended in Grynkiv and Stentoft (2018).

Likewise, assume the variance of the error is constant. If we consider a mean corrected Poisson so that $\psi \left(u\right)=\exp \left(\theta \left(\exp \left(iu\right)-1-iu\right)\right)$ with mean parameter $\theta$, we can identify the distribution of $p_t$ as an infinite weighted sum of Poisson random variables of increasing means (j + 1) $\theta$, but whose relative importance declines with a power of $\pi .$

This case can be extended to include intercepts, in which case,

$$\begin{gathered} \varphi \left(u\right)=(\exp \left(i{\alpha }_{1}u\right)\varphi \left(u\right)\pi +\left(\exp \left(i{\alpha }_{2}u\right)\left(1-\pi \right)\right)\psi \left(u\right) \\ \varphi \left(u\right)=\frac{\exp \left(i{\alpha }_{2}u\right)\left(1-\pi \right)\psi \left(u\right)}{1-\exp \left(i{\alpha }_{1}u\right)\pi \psi \left(u\right)} \end{gathered}$$

Since $\left|\exp \left(i{\alpha }_{1}u\right)\pi \psi \left(u\right)\right|<1,$ the analysis proceeds as before, and:

$$\varphi \left(u\right)=\left(1-\pi \right)\exp \left(i{\alpha }_{2}u\right){\displaystyle \sum }_{j=0}^{\infty }{\pi }^j\psi {\left(u\right)}^{j+1}\exp \left(ij{\alpha }_{1}u\right)$$

We see that the jth component is as above, but has a mean augmented by $j{\alpha }_1+{\alpha }_2$. For other examples, we refer the interested reader to Pourahmadi (1988) Section 4, who derived marginal distributions for different processes.

The results can be generalized to Vector threshold auto-regressions. These results are not presented here but are available upon request. The interested reader may also refer to Grynkiv and Stentoft (2018) for some discussion.

3. A Caveat on Simulating the Steady State

Simulating the steady state for processes similar to those considered in Section 2 is not as straight forward as it may appear at first glance and warrants further consideration, which this section seeks to provide. Generating a discrete Markov chain, which is essentially a variable that takes discrete integer values of 0 or 1 depending on the transition matrix P, does not generate a steady-state Markov chain, but rather a path to the steady state. It is common in this literature to simulate steady-state paths rather than a steady-state Markov chain. Indeed, in our earlier work, we worked with paths and not steady states, as did Timmerman in his article (see Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6 in his articles for example). While simulating steady-state paths sufficed for our earlier work, we need to simulate the steady state in order to corroborate our results from Section 2. Otherwise, the underlying moments of the simulated series can be different even if they converge to the same steady state.

This path is dependent upon the transition matrix P. If states are persistent, as determined by the transition probabilities of the process staying in the prevailing state ($p_{ii})$, this path may diverge significantly from the steady state. On the opposite spectrum is a transition matrix with frequent state switches, due to higher switching probabilities, which will take a different path to the steady state. Although the steady-state probabilities of both paths are identical, the dynamics vary due to the different paths taken by the processes.

We need to consider how discrete Markov chains converge to their steady states. The usual definition is based on total variation distance, and considers the supremum, taken over measurable subsets A, of the absolute difference between $v\left(A\right)$ and $u\left(A\right)$, where $v()$ and $u()$ are the two probability measures (see Rosenthal 1995). Whilst our process will converge in this sense, it will almost surely not converge along a sample path. Intuitively, it keeps moving from state to state.

We illustrate this by considering two different transition matrices that correspond to the steady-state probability vector $\left[\begin{array}{c}0.5\\ 0.5\end{array}\right]$. We consider a transition matrix with highly persistent states $P_1=\left(\begin{array}{cc}0.9& 0.1\\ 0.1& 0.9\end{array}\right)$ and a transition matrix with frequent state switches $P_2=\left(\begin{array}{cc}0.5& 0.5\\ 0.5& 0.5\end{array}\right)$. While both processes approach the same steady states, the simulated series have different probability distributions. Specifically, persistence of the non-stationary process, corresponding to ${\beta }_j>1$, causes the path to diverge far from the steady-state values, resulting in a process that has extreme values with significant probability, which also obtains very high kurtosis. With frequent state switches, the simulated series comes closer to a normally distributed process.

To remedy this, each of the above paths was simulated for 10,000 periods and the process was repeated 2000 times. The parameter values in the two switching states are 0 in state 1 and 1 in state 2. We also carried out simulations when the switching states corresponded to values of 0.1 in state 1 and 1.1 in state 2. We recorded the average of the first four moments of both paths in

Table 1. Moments of TAR(1) with Markov-switching exogenous variables.

Steady State Vector	Transition Matrix	$\mathit{\delta }$	Mean	Stdev	Skewness	Kurtosis
[0.50 0.50]	$P_1=\left(\begin{array}{cc}0.9& 0.1\\ 0.1& 0.9\end{array}\right)$	0	0.0004	2.450	0.00289	8.786
[0.50 0.50]	$P_1=\left(\begin{array}{cc}0.9& 0.1\\ 0.1& 0.9\end{array}\right)$	0.10	−0.170	78.17	0.2761	424.7
[0.50 0.50]	$P_2=\left(\begin{array}{cc}0.5& 0.5\\ 0.5& 0.5\end{array}\right)$	0	−0.0001	1.4139	−0.0004	4.480
[0.50 0.50]	$P_2=\left(\begin{array}{cc}0.5& 0.5\\ 0.5& 0.5\end{array}\right)$	0.10	0.0000	1.6014	−0.0012	6.837

Table 1 reports the first four moments of simulated threshold autoregressive series with Markov-switching exogenous variables. Each series is 10,000 observations long, and 2000 series were simulated for each set of values. P represents the transition matrix and delta represents a parameter that determines the values of ${\beta }_j$ in the two states; ${\beta }_1=\delta ,{\beta }_2=1+\delta$.

. As mentioned previously, the persistent path obtains a much higher kurtosis and standard deviation than the more frequently switching path, even though both paths continue to be symmetrically distributed. An alternative approach would be to write out the solution to Equation (1) and simulate directly by taking long samples of the error term and the exogenous process. From Table 1, we observe that the highest 2nd and 4th moments are obtained for the persistent state when ${\beta }_1=0.1 \mathrm{and} {\beta }_2=1.1$. This is despite the fact that this process also converges to a steady-state vector $\left[0.50 0.50\right]$. The moments are significantly different for the alternative paths corresponding to $P_2$.

Figure 1 and Figure 2 below highlight the different distributions that result from the different paths, along with cumulative probabilities corresponding to the normal probability density function’s quintile values. Tail probabilities are much higher for the more persistent path, which also has more extreme values, especially when the non-stationary state becomes explosive. Tail probabilities are 28.7% and 30.8% for $\delta =0$ and $\delta =0.10$, respectively, in Figure 1, which correspond to the more persistent transition matrix, $P_1$. In fact, the distribution corresponding to $\delta =0.10$ appears like a horizontal line instead of a bell-shaped curve.

On the other hand, the distributions of simulated series corresponding to the more frequently switching path (transition matrix $P_2$) have lower tail probabilities and distributions that are closer to the standard normal distribution. Figure 2 contrasts the distribution for the series corresponding to path $P_2$ to a standard normal. While this simulated distribution has heavier tails, as evident from the higher probabilities corresponding to normal quintiles, its 2nd and 4th moments are much closer to the normal than to the simulated series for path $P_1$. Similarly, for the distributions in Figure 2, corresponding to ${\beta }_1=0.10 and {\beta }_2=1.10$, path $P_2$ look much closer to a normal distribution than to their $P_1$ counterparts in Figure 1.

Thus, we need a different approach to generate the steady-state distribution of the threshold autoregressive process with Markovian triggers that are independent of the transition matrix, subject to the same steady state. It is important to understand that the results derived in the section above correspond to the steady state itself and not to the path of the process tending to a steady state, which, as we have shown in this section, depends on the transition matrix. We describe how we simulate the steady-state distribution in the next section.

4. Simulation Results

In order to simulate a distribution for the steady state, we observed the Markov chain 10,000 times. However, each observation was 1000 time periods or steps apart. Thus, the Markov chain we simulated was 10,000,000 periods long, and the steady-state simulation was 10,000 observations in length. We verified that each series simulated this way converges to the steady-state probability vector while at the same time being independent of the transition matrix, P. We did this by simulating steady states for different transition matrices P that shared the same steady states.1 Thus, the analysis considered in this section depends only on the steady-state vector and not on the transition matrices. The steady-state chain was then used to simulate the following threshold model:

$$y_t=\alpha +{\beta }_{it}{y}_{t-1}+{\eta }_t$$

where ${\beta }_{it}$ depends on the value taken by the exogenous discrete Markov state variable $z_i$. Our simulations assumed that ${\eta }_t~N\left(0,1\right)$ and that $\alpha =0$. In order for our simulated series to have a steady-state stationary distribution, we required that the criterion ${{\displaystyle \sum }}_{i=1}^2{\pi }_i\ln \left|{\beta }_i\right|<0$ be satisfied. Since we considered a two-state process, the criterion can alternatively be written as ${\beta }_1^{\pi }{\beta }_2^{1-\pi }<1$. The criterion is trivially satisfied when ${\beta }_1=0$ as $\ln \left(0\right)=-\infty$. When ${\beta }_1=\delta$ and ${\beta }_2=\left(1+\delta \right)$, the criterion becomes ${\delta }^{\pi }{\left(1+\delta \right)}^{1-\pi }<1$. Note, that the expression ${\delta }^{\pi }{\left(1+\delta \right)}^{1-\pi }$ is maximized for $\pi =0$.

We considered a maximum $\delta$ of 0.10, and checked that the criterion was satisfied for all our simulations and that a steady-state distribution does exist. For exposition, we included the value taken by the criterion function for each set of simulations in column 7 of

Table 2. Average moments of TAR(1) with changing Steady State Vectors.

Steady State Vector $\left[\mathit{\pi } 1-\mathit{\pi }\right]$	$\mathit{\delta }$	Mean	Stdev	Skewness	Kurtosis	$\mathbf{Criteria}={{\displaystyle \sum }}_{\mathit{i}=1}^2{\mathit{\pi }}_{\mathit{i}}\mathbf{ln}\left|{\mathit{\beta }}_{\mathit{i}}\right|$
[0.90 0.10]	0	0.0001	1.0544	−0.0012	3.306	$-\infty$
[0.90 0.10]	0.01	−0.0007	1.0552	0.0004	3.311	−4.1437
[0.90 0.10]	0.03	0.0001	1.0583	−0.0008	3.336	−3.1529
[0.90 0.10]	0.05	0.0007	1.0608	−0.0006	3.362	−2.6913
[0.90 0.10]	0.10	0.0007	1.0714	−0.0007	3.442	−2.0628
[0.80 0.20]	0	0.0004	1.118	−0.002	3.597	$-\infty$
[0.80 0.20]	0.01	−0.0000	1.121	−0.0007	3.630	−3.6821
[0.80 0.20]	0.03	0.0001	1.128	−0.0016	3.705	−2.7993
[0.80 0.20]	0.05	0.0005	1.134	0.003	3.762	−2.3868
[0.80 0.20]	0.10	0.0001	1.152	0.0003	3.967	−1.8230
[0.60 0.40]	0	0.0004	1.291	−0.003	4.195	$-\infty$
[0.60 0.40]	0.01	−0.0011	1.300	−0.001	4.277	−2.7591
[0.60 0.40]	0.03	−0.0001	1.316	−0.005	4.421	−2.0921
[0.60 0.40]	0.05	−0.0007	1.339	−0.002	4.660	−1.7779
[0.60 0.40]	0.10	−0.0004	1.400	0.002	5.410	−1.3434
[0.50 0.50]	0	0.0017	1.413	0.0026	4.470	$-\infty$
[0.50 0.50]	0.01	−0.0004	1.428	−0.0042	4.612	−2.2976
[0.50 0.50]	0.03	−0.0006	1.460	−0.0001	4.915	−1.7385
[0.50 0.50]	0.05	−0.0004	1.494	−0.0021	5.302	−1.4735
[0.50 0.50]	0.10	0.001	1.601	0.0227	6.720	−1.1036

Table 2 reports the first four moments of the simulated series along with the stationary criterion. Each row corresponds to 2000 threshold autoregressive simulations with Markov triggers, 10,000 observations long; columns 3–6 report average moments.

below.

Steady-state distributions obtained this way can be analyzed through the results derived in Section 2. For each set of steady-state vectors in Table 2, we simulated threshold autoregressive series (as described above) 2000 times. The parameter values in the two switching states are $\delta$ in state 1 and $1+\delta$ in state 2. Thus, the first state is always stationary, while the second state is either a random walk or explosive.

Some patterns are exhibited clearly. The distributions of the series appear to be centered on zero, and statistically their skewness (column 5) is not significantly different from zero, which follows from Theorem 3. Since the disturbance term has an even distribution, it follows that the distribution of our simulated series will also be even, symmetric, and centered on zero, i.e., $\psi \left(u\right)=\psi \left(-u\right)$ implies that $\varphi \left(u\right)=\varphi \left(-u\right)$.

The standard deviation (column 2), however, does appear to be much larger than the standard deviation of the error process driving the threshold process. The series have excess kurtosis (column 6), which should not come as a surprise since the series display non-stationary behavior when ${\beta }_i\ge 1$. This state leads to excess kurtosis and higher standard deviation. As we deviate from our base case, ($\delta =0)$, we note a clear pattern in the 2nd and 4th moments of the series. Both the standard deviation and kurtosis start to increase, since this behavior is caused by the non-stationary state becoming increasingly more explosive (${\beta }_i=1+\delta , when \delta >0)$. The pattern is repeated irrespective of the steady-state vector chosen.

Unsurprisingly, the higher the steady-state probability of the stationary state (${\beta }_i=\delta$), the closer the process’s distribution is to a normal distribution. This is reflected in the first four moments. For instance, when the stationary state occurs 90% of the time (as in rows 2–6 of Table 2), the standard deviation and kurtosis are both lower compared to corresponding cases (i.e., same $\delta$) when the stationary state occurs less than 90% of the time. When $\delta =0.05$, the standard deviation and kurtosis respectively are 1.061 and 3.362 for $\pi =0.90$, 1.134 and 3.762 for $\pi =0.80$, 1.339 and 4.660 for $\pi =0.60$, and 1.494 and 5.302 for $\pi =0.5$. While all distributions appear symmetric and centered on 0, they become increasingly leptokurtic as $\pi$ falls.

We also plotted sample distributions for $\delta =0$ and $\delta =0.10$ for each set of steady-state vectors (Figure 3, Figure 4, Figure 5 and Figure 6) and carried out quintile analysis by calculating the weights in the distribution corresponding to the quintile values of a normal distribution, i.e., we found $P\left(y_t\right)<q_1, q_1<P\left(y_t\right)<q_2,q_2<P\left(y_t\right)<q_3, q_3<P\left(y_t\right)<q_4,\mathrm{and}P\left(y_t\right)>q_4$, where $q_i$ corresponds to the normal distribution’s quintile values. We note that tail probabilities, i.e., those corresponding to the 1st and 5th normal quintile values, go up as the steady-state probabilities for the non-stationary state go up. They are also dependent on the value of $\delta$, and as we increase explosivity in the non-stationary state (by increasing $\delta$), tail probabilities increase and the distributions moves farther away from a normal. The symmetry of the distribution is also reflected in these probabilities.

The results above depend only on the steady-state vector $\left[\begin{array}{c}\pi \\ 1-\pi \end{array}\right]$ and not on the transition matrix (the dynamic path) that generates this steady-state vector. We note that all series generated this way have excess kurtosis due to the presence of a non-stationary state.2 Below, we simulate probability density functions for some of the cases considered in Table 2. The graphs only report the distribution for $\delta =0$ and $\delta =0.10$, so we can analyze how far the distribution moves from a normal as we increase $\delta$ and increase the probability of the non-stationary part of the distribution. We also draw a comparison with a normal distribution for illustrative purposes.

It is worth mentioning at this point that $\delta =0.10$ may not present as interesting an empirical case as say $\delta =0.01$ or $\delta =0.03$. In empirical work (e.g., Ahmed and Satchell 2018), estimates have suggested that even with an explosive state in a threshold autoregressive model, the $\delta$ of the explosive state from the random walk state does not extend beyond $0.03$. Our analysis is thus meant to provide a graphical view of how the distribution evolves as $\delta$ is varied, and this is most obvious when considering the two extreme values. Distributional analysis for other values was also considered and is available upon request. While results for other values of $\delta$ follow in the same direction as those considered here, the differences are less stark as would be expected. For researchers interested in modelling non-stationary behavior using threshold auto-regressions, a value of $\delta <0.03$ may be more appropriate.

From Figure 3a,b, we can see that when the stationary state is dominant $\left(\pi =0.90\right)$, the distribution appears very close to a normal distribution, even if $\delta$ is increased from 0 to 0.10. Contrast this distribution with Figure 6a,b, where the stationary and non-stationary states occur with equal probability. The distributions in Figure 6a,b are significantly leptokurtic, with tail probabilities of 25.2% and 26.2%, respectively, as opposed to tail probabilities of 20.8% and 21% in Figure 3a,b, respectively.

These results are interesting, particularly for those relying on DSGE models. In finance, asset prices often exhibit locally non-stationary behavior, which leads to leptokurtosis in the series. Similarly, macroeconomists often consider different shock mechanisms in DSGE models. The results in this article will assist macroeconomists in considering shock processes that follow threshold auto-regressions with a non-stationary state. Since we have outlined a procedure for deriving analytical expressions, this would enable macroeconomists to analyze locally explosive shock processes which nevertheless are stationary overall and facilitate the implementation of numerical solutions.

5. Conclusions

In this article, we have derived formulae for characteristic functions for threshold autoregressive models of order 1, which have a Markovian state-switching variable. In doing so, we have not only improved the results first considered in Timmermann (2000), but have also generalized the formulae to a great degree. These formulae will allow readers, if they are so inclined, to derive analytical moments for TAR models in this class for a range of different error specifications. We believe that these will have applications in both finance and applied macroeconomics, given increasing interest in these models.

Considering a special case of interest, we have also shown using simulation analysis that the existence of a non-stationary state in a TAR model can cause the distribution of a particular series to deviate significantly from normality. The further away the non-stationary state moves from a random walk, the farther away the distribution is from that of a normal. Models for asset prices often consider a mixture of stationary and non-stationary states, and we believe that this simulation analysis will aid researchers in better understanding the behavior of asset prices that go through locally explosive states. In Section 1 and Section 2, we have further outlined how our results may be applicable to applied macroeconomic and finance literatures.

Indeed, articles like Ahmed and Satchell (2018) would have found these results useful. While Ahmed and Satchell (2018) provided an empirical example to estimate market efficiency using a threshold autoregressive model, they were unable to calculate means and variances for their series. Through the results in our article, they would have been able to compare the two asset price series using well-known metrics in finance, such as Sharpe ratios. Similarly, any researchers modelling asset prices through threshold autoregressive models with a Markov-switching variable will be able to calculate moments and provide a more holistic narrative.

Furthermore, the contributions in the article will allow researchers to consider non-normal distributions for asset price series. Normality or log-normality is often assumed in the asset pricing literature to allow derivation of analytical expressions at the cost of accuracy. As we have shown, not only are threshold autoregressive models better at capturing asset pricing-like behavior, but they also allow for analytical solutions that can be useful in helping us better understand asset pricing behavior.

Another avenue worth considering from a statistical point of view will be to consider relevant moment-generating functions for the error process, and analyzing the resulting distributions for the underlying process, which we have referred to as the price process. Whilst we have considered the standard normal and Poisson distributions, future research could consider other distributions, which may inform future Finance and Economics research.

Thus, our article makes significant contributions to the existing literature on TAR models and offers insights into how such models may be used in applied economics and finance.