A mixture autoregressive model based on Student's $t$-distribution

A new mixture autoregressive model based on Student's $t$-distribution is proposed. A key feature of our model is that the conditional $t$-distributions of the component models are based on autoregressions that have multivariate $t$-distributions as their (low-dimensional) stationary distributions. That autoregressions with such stationary distributions exist is not immediate. Our formulation implies that the conditional mean of each component model is a linear function of past observations and the conditional variance is also time varying. Compared to previous mixture autoregressive models our model may therefore be useful in applications where the data exhibits rather strong conditional heteroskedasticity. Our formulation also has the theoretical advantage that conditions for stationarity and ergodicity are always met and these properties are much more straightforward to establish than is common in nonlinear autoregressive models. An empirical example employing a realized kernel series based on S&P 500 high-frequency data shows that the proposed model performs well in volatility forecasting.


Introduction
Different types of mixture models are in widespread use in various fields. Overviews of mixture models can be found, for example, in the monographs of McLachlan & Peel (2000) and Frühwirth-Schnatter (2006). In this paper, we are concerned with mixture autoregressive models that were introduced by Le et al. (1996) and further developed by Wong & Li (2000, 2001a (for further references, see ).
In mixture autoregressive models the conditional distribution of the present observation given the past is a mixture distribution where the component distributions are obtained from linear autoregressive models. The specification of a mixture autoregressive model typically requires two choices: choosing a conditional distribution for the component models and choosing a functional form for the mixing weights. In a majority of existing models a Gaussian distribution is assumed whereas, in addition to constants, several different time-varying mixing weights (functions of past observations) have been considered in the literature.
Instead of a Gaussian distribution,  proposed using Student's t-distribution. A major motivation for this comes from the heavier tails of the t-distribution which allow the resulting model to better accommodate for the fat tails encountered in many observed time series, especially in economics and finance. In the model suggested by , the conditional mean and conditional variance of each component model are the same as in the Gaussian case (a linear function of past observations and a constant, respectively), and what changes is the distribution of the independent and identically distributed error term: instead of a standard normal distribution, a Student's t-distribution is used. This is a natural approach to formulate the component models and hence also a mixture autoregressive model based on the t-distribution.
In this paper, we also consider a mixture autoregressive model based on Student's t-distribution, but our specification differs from that used by . Our starting point is the characteristic feature of linear Gaussian autoregressions that stationary distributions (of consecutive observations) as well as conditional distributions are Gaussian. We imitate this feature by using a (multivariate) Student's t-distribution and, as a first step, construct a linear autoregression in which both conditional and (low-dimensional) stationary distributions have Student's t-distributions. This leads to a model where the conditional mean is as in the Gaussian case (a linear function of past observations) whereas the conditional variance is no longer constant but depends on a quadratic form of past observations. These linear models are then used as component models in our new mixture autoregressive model which we call the StMAR model.
Our StMAR model has some very attractive features. Like the model of , it can be useful for modelling time series with regime switching, multimodality, and conditional heteroskedasticity. As the conditional variances of the component models are time-varying, the StMAR model can potentially accommodate for stronger forms of conditional heteroskedasticity than the model of . Our formulation also has the theoretical advantage that, for a pth order model, the stationary distribution of p + 1 consecutive observations is fully known and is a mixture of particular Student's t-distributions. Moreover, stationarity and ergodicity are simple consequences of the definition of the model and do not require complicated proofs.
Finally, a few notational conventions. All vectors are treated as column vectors and we write x = (x 1 , . . . , x n ) for the vector x where the components x i may be either scalars or vectors. The notation X ∼ n d (µ, Γ) signifies that the random vector X has a d-dimensional Gaussian distribution with mean µ and (positive definite) covariance matrix Γ. Similarly, by X ∼ t d (µ, Γ, ν) we mean that X has a d-dimensional Student's t-distribution with mean µ, (positive definite) covariance matrix Γ, and degrees of freedom ν (assumed to satisfy ν > 2); the density function and some properties of the multivariate Student's t-distribution employed are given in an Appendix. The notation 1 d is used for a d-dimensional vector of ones, ı d signifies the vector (1, 0, . . . , 0) of dimension d, and the identity matrix of dimension d is denoted by I d . The Kronecker product is denoted by ⊗, and vec(A) stacks the columns of matrix A on top of one another.

Linear Student's t autoregressions
In this section we briefly consider linear pth order autoregressions that have multivariate Student's t-distributions as their stationary distributions. First, for motivation and to develop notation, consider a linear Gaussian autoregression z t (t = 1, 2, . . .) generated by where the error terms e t are independent and identically distributed with a standard normal distribution, and the parameters satisfy ϕ 0 ∈ R, ϕ = (ϕ 1 , . . . , ϕ p ) ∈ S p , and σ > 0, where is the stationarity region of a linear pth order autoregression. Denoting z t = (z t , . . . , z t−p+1 ) and z + t = (z t , z t−1 ), it is well known that the stationary solution z t to (1) satisfies where the last relation defines the conditional distribution of z t given z t−1 and the quantities Γ p , γ 0 , γ p , µ, and Γ p+1 are defined via Two essential properties of linear Gaussian autoregressions are that they have the distributional features in (3) and the representation in (1). It is not immediately obvious that linear autoregressions based on Student's t-distribution with similar properties exist (such models have, however, appeared at least in Spanos (1994) and Heracleous & Spanos (2006)). Suppose that for a random vector in R p+1 it holds that (z, z) ∼ t p+1 (µ1 p+1 , Γ p+1 , ν) where ν > 2 (and other notation is as above in (4)). Then (for details, see the Appendix) the conditional distribution of z given z is z | z ∼ t 1 (µ(z), σ 2 (z), ν + p), where We now state the following theorem (proofs of all theorems are in the Supplementary Material).
Results (i) and (ii) in Theorem 1 are comparable to properties (3) and (1) in the Gaussian case. Part (i) shows that both the stationary and conditional distributions of z t are t-distributions, whereas part (ii) clarifies the connection to standard AR(p) models. In contrast to linear Gaussian autoregressions, in this t-distributed case z t is conditionally heteroskedastic and has an 'AR(p)-ARCH(p)' representation (here ARCH refers to autoregressive conditional heteroskedasticity).

A mixture autoregressive model based on Student's t-distribution
3.1 Mixture autoregressive models Let y t (t = 1, 2, . . .) be the real-valued time series of interest, and let F t−1 denote the σ-algebra generated by {y t−j , j > 0}. We consider mixture autoregressive models for which the conditional density function of y t given its past, f (· | F t−1 ), is of the form where the (positive) mixing weights α m,t are F t−1 -measurable and satisfy M m=1 α m,t = 1 (for all t), and the f m (· | F t−1 ), m = 1, . . . , M , describe the conditional densities of M autoregressive component models. Different mixture models are obtained with different specifications of the mixing weights α m,t and the conditional densities f m (· | F t−1 ).
Starting with the specification of the conditional densities f m (· | F t−1 ), a common choice has been to assume the component models to be linear Gaussian autoregressions. For the mth component model (m = 1, . . . , M ), denote the parameters of a pth order linear autoregression with ϕ m,0 ∈ R, ϕ m = (ϕ m,1 , . . . , ϕ m,p ) ∈ S p , and σ m > 0. Also set y t−1 = (y t−1 , . . . , y t−p ). In the Gaussian case, the conditional densities in (8) take the form (m = 1, . . . , M ) where φ(·) signifies the density function of a standard normal random variable, µ m,t = ϕ m,0 +ϕ ′ m y t−1 is the conditional mean function (of component m), and σ 2 m > 0 is the conditional variance (of component m), often assumed to be constant. Instead of a Gaussian density,  consider the case where f m (· | F t−1 ) is the density of Student's t-distribution with conditional mean and variance as above, µ m,t = ϕ m,0 + ϕ ′ m y t−1 and a constant σ 2 m , respectively. In this paper, we also consider a mixture autoregressive model based on Student's t-distribution, but our formulation differs from that used by . In Theorem 1 it was seen that linear autoregressions based on Student's t-distribution naturally lead to the conditional distribution t 1 (µ(·), σ 2 (·), ν + p) in (6). Motivated by this, we consider a mixture autoregressive model in which the conditional densities f m (y t | F t−1 ) in (8) are specified as where the expressions for µ m,t = µ m (y t−1 ) and σ 2 m,t = σ 2 m (y t−1 ) are as in (5) except that z is replaced with y t−1 and all the quantities therein are defined using the regime specific parameters ϕ m,0 , ϕ m , σ m , and ν m (whenever appropriate a subscript m is added to previously defined notation, e.g., µ m or Γ m,p ).
In this paper, we propose mixing weights that are similar to those used by Glasbey (2001) and . Specifically, we set where the α m ∈ (0, 1), m = 1, . . . , M , are unknown parameters satisfying M m=1 α m = 1. Note that the Student's t density appearing in (11) corresponds to the stationary distribution in Theorem 1(i): If the y t 's were generated by a linear Student's t autoregression described in Section 2 (with a subscript m added to all the notation therein), the stationary distribution of y t−1 would be characterized by t p (y t−1 ; µ m 1 p , Γ m,p , ν m ). Our definition of the mixing weights in (11) is different from that used in Glasbey (2001) and  in that these authors employed the n p (y t−1 ; µ m 1 p , Γ m,p ) density (corresponding to the stationary distribution of a linear Gaussian autoregression) instead of the Student's t density t p (y t−1 ; µ m 1 p , Γ m,p , ν m ) we use.

The Student's t mixture autoregressive model
Equations (8), (9), and (11) define a model we call the Student's t mixture autoregressive, or StMAR, model. When the autoregressive order p or the number of mixture components M need to be emphasized we refer to an StMAR(p,M ) model. We collect the unknown parameters of an StMAR model in the vector θ = (ϑ 1 , . . . , ϑ M , α 1 , . . . , α M −1 ) ((M (p + 4) − 1) × 1), where ϑ m = (ϕ m,0 , ϕ m , σ 2 m , ν m ) (with ϕ m ∈ S p , σ 2 m > 0, and ν m > 2) contains the parameters of each component model (m = 1, . . . , M ) and the α m 's are the parameters appearing in the mixing weights (11); the parameter α M is not included due to the restriction M m=1 α m = 1. The StMAR model can also be presented in an alternative (but equivalent) form. To this end, let P t−1 (·) signify the conditional probability of the indicated event given F t−1 , and let ε m,t be a sequence of independent and identically distributed random variables with a t 1 (0, 1, ν m + p) distribution such that ε m,t is independent of {y t−j , j > 0} (m = 1, . . . , M ). Furthermore, let s t = (s 1,t , . . . , s M,t ) be a sequence of (unobserved) M -dimensional random vectors such that, conditional on F t−1 , s t and ε m,t are independent (for all m). The components of s t are such that, for each t, exactly one of them takes the value one and others are equal to zero, with conditional probabilities P t−1 (s m,t = 1) = α m,t , m = 1, . . . , M . Now y t can be expressed as where σ m,t is as in (9). This formulation suggests that the mixing weights α m,t can be thought of as (conditional) probabilities that determine which one of the M autoregressive components of the mixture generates the observation y t . It turns out that the StMAR model has some very attractive theoretical properties; the carefully chosen conditional densities in (9) and the mixing weights in (11) are crucial in obtaining these properties. The following theorem shows that there exists a choice of initial values y 0 such that y t is a stationary and ergodic Markov chain. Importantly, an explicit expression for the stationary distribution is also provided.
The stationary distribution of y t is a mixture of M p-dimensional t-distributions with constant mixing weights α m . Hence, moments of the stationary distribution of order smaller than min (ν 1 , . . . , ν M ) exist and are finite. As can be seen from the proof of Theorem 2 (in the Supplementary Material), the stationary distribution of the vector (y t , y t−1 ) is also a mixture of M t-distributions with density of the same form, M m=1 α m t p+1 (µ m 1 p+1 , Γ m,p+1 , ν m ). Thus the mean, variance, and first p autocovariances of y t are (here the connection between γ m,j and Γ m,p+1 is as in (4)) Subvectors of (y t , y t−1 ) also have stationary distributions that belong to the same family (but this does not hold for higher dimensional vectors such as (y t+1 , y t , y t−1 )).
The fact that an explicit expression for the stationary (marginal) distribution of the StMAR model is available is not only convenient but also quite exceptional among mixture autoregressive models or other related nonlinear autoregressive models (such as threshold or smooth transition models). Previously, similar results have been obtained by Glasbey (2001) and  in the context of mixture autoregressive models that are of the same form but based on the Gaussian distribution (for a few rather simple first order examples involving other models, see Tong (2011, Section 4.2)).
From the definition of the model, the conditional mean and variance of y t are obtained as Except for the different definition of the mixing weights, the conditional mean is as in the Gaussian mixture autoregressive model of . This is due to the well-known fact that in the multivariate t-distribution the conditional mean is of the same linear form as in the multivariate Gaussian distribution. However, unlike in the Gaussian case, the conditional variance of the multivariate t-distribution is not constant. Therefore, in (13) we have the time-varying variance component σ 2 m,t which in the models of  and  is constant (in the latter model the mixing weights are also constants). In (13) both the mixing weights α m,t and the variance components σ 2 m,t are functions of y t−1 , implying that the conditional variance exhibits nonlinear autoregressive conditional heteroskedasticity. Compared to the aforementioned previous models our model may therefore be useful in applications where the data exhibits rather strong conditional heteroskedasticity.

Estimation
The parameters of an StMAR model can be estimated by the method of maximum likelihood (details of the numerical optimization methods employed and of simulation experiments are available in the Supplementary Material). As the stationary distribution of the StMAR process is known it is even possible to make use of initial values and construct the exact likelihood function and obtain exact maximum likelihood estimates. Assuming the observed data y −p+1 , . . . , y 0 , y 1 , . . . , y T and stationary initial values, the log-likelihood function takes the form where An explicit expression for the density appearing in (15) is given in (10), and the notation for µ m,t and σ 2 m,t is explained after (9). Although not made explicit, α m,t , µ m,t , and σ 2 m,t , as well as the quantities µ m , γ m,p , and Γ m,p , depend on the parameter vector θ.
In (14) it has been assumed that the initial values y 0 are generated by the stationary distribution. If this assumption seems inappropriate one can condition on initial values and drop the first term on the right hand side of (14). In what follows we assume that estimation is based on this conditional log-likelihood, namely L (c) which we, for convenience, have also scaled with the sample size. Maximizing L (c) T (θ) with respect to θ yields the maximum likelihood estimator denoted byθ T .
The permissible parameter space of θ, denoted by Θ, needs to be constrained in various ways. The stationarity conditions ϕ m ∈ S p , the positivity of the variances σ 2 m , and the conditions ν m > 2 ensuring existence of second moments are all assumed to hold (for m = 1, . . . , M ). Throughout we assume that the number of mixture components M is known, and this also entails the requirement that the parameters α m (m = 1, . . . , M ) are strictly positive (and strictly less than unity whenever M > 1). Further restrictions are required to ensure identification. Denoting the true parameter value by θ 0 and assuming stationary initial values, the condition needed is that l t (θ) = l t (θ 0 ) almost surely only if θ = θ 0 . An additional assumption needed for this is From a practical point of view this assumption is not restrictive because what it essentially requires is that the M component models cannot be 'relabeled' and the same StMAR model obtained. We summarize the restrictions imposed on the parameter space as follows. Asymptotic properties of the maximum likelihood estimator can now be established under conventional high-level conditions. Denote Theorem 3. Suppose y t is generated by the stationary and ergodic StMAR process of Theorem 2 and that Assumption 1 holds. Thenθ T is strongly consistent, i.e.,θ T → θ 0 almost surely. Suppose further that Of the conditions in this theorem, (i) states that a central limit theorem holds for the score vector (evaluated at θ 0 ) and that the information matrix is positive definite, (ii) is the information matrix equality, and (iii) ensures the uniform convergence of the Hessian matrix (in some neighbourhood of θ 0 ). These conditions are standard but their verification may be tedious.
Theorem 3 shows that the conventional limiting distribution applies to the maximum likelihood estimatorθ T which implies the applicability of standard likelihood-based tests. It is worth noting, however, that here a correct specification of the number of autoregressive components M is required. In particular, if the number of component models is chosen too large then some parameters of the model are not identified and, consequently, the result of Theorem 3 and the validity of the related tests break down. This particularly happens when one tests for the number of component models. Such tests for mixture autoregressive models with Gaussian conditional densities (see (8)) are developed by Meitz & Saikkonen (2017). The testing problem is highly nonstandard and extending their results to the present case is beyond the scope of this paper.
Instead of formal tests, in our empirical application we use information criteria to infer which model fits the data best. Similar approaches have also been used by  and others. Note that once the number of regimes is (correctly) chosen, standard likelihood-based inference can be used to choose regime-wise autoregressive orders and to test other hypotheses of interest.

Empirical example
Modeling and forecasting financial market volatility is key to manage risk. In this application we use the realized kernel of Barndorff-Nielsen et al. (2008) as a proxy for latent volatility. We obtained daily realized kernel data over the period 3 January 2000 through 20 May 2016 for the S&P 500 index from the Oxford-Man Institute's Realized Library v0.2 (Heber et al., 2009). Figure 1 shows the in-sample period (Jan 3, 2000-June 3, 2014; 3597 observations) for the S&P 500 realized kernel data (RK t ), which is nonnegative with a distribution exhibiting substantial skewness and excess kurtosis (sample  Table 1 (dot-dash) for the log(RK t ) series. The mixing weightsα 1,t are scaled from (0, 1) to (min log(RK t ), max log(RK t )). Right panel: A kernel density estimate of the log(RK t ) observations (solid), and the mixture density (dashes) implied by the same StMAR model as in the left panel. skewness 14.3, sample kurtosis 380.8). We follow the related literature which frequently use logarithmic realized kernel (log(RK t )), to avoid imposing additional parameter constraints, and to obtain a more symmetric distribution, often taken to be approximately Gaussian. The log(RK t ) data, also shown in Figure 1, has a sample skewness of 0.5 and kurtosis of 3.5. Visual inspection of the time series plots of the RK t and log(RK t ) data suggests that the two series exhibit changes at least in levels and potentially also in variability. A kernel estimate of the density function of the log(RK t ) series also suggest the potential presence of multiple regimes. Table 1 reports estimation results for three selected StMAR models (for further details, see the Supplementary Material). Following Wong & Li (2001a), Li et al. (2015), we use information criteria for model comparison. For the log(RK t ) data in-sample period the Akaike information criterion (aic) favours the StMAR(4,3) model, the Hannan-Quinn information criterion (hqc) the StMAR(4,2) model, and the Bayesian information criterion (bic) the simpler StMAR(4,1) model. In view of the approximate standard errors in Table 1, the estimation accuracy appears quite reasonable except for the degrees of freedom parameters. Taking the sum of the autoregressive parameters as a measure of persistence, we find that the estimated persistence for the first regime of the StMAR(4,2) is 0.909 and 0.489 for the second regime, suggesting that persistence is rather strong in the first regime and moderate in the second regime.
Numerous alternative models for volatility proxies have been proposed. We employ Corsi's (2009) heterogeneous autoregressive (HAR) model as it is arguably the most popular reference model for forecasting proxies such as the realized kernel. We also consider a pth-order autoregression as the AR(p) often performs well in volatility proxy forecasting. The StMAR models are estimated using maximum likelihood, and the reference AR and HAR models by ordinary least squares. We use a fixed scheme, where the parameters of our volatility models are estimated just once using data from Jan 3, 2000-June 3, 2014. These estimates are then used to generate all forecasts. The remaining 496 observations of our sample are used to compare the forecasts from the alternative models. As discussed in , computing multi-step-ahead forecasts for mixture models like the StMAR is rather complicated. For this reason we use computer driven forecasts to predict future volatility: For each out-of-sample date T , and for each alternative model, we simulate 500,000 sample paths. Each  path is of length 22 (representing one trading month) and conditional on the information available at date T . In these simulations unknown parameters are replaced by their estimates. As the simulated paths are for log(RK t ), and our object of interest is RK t , an exponential transformation is applied. We examine daily, weekly (5 day), biweekly (10 day), and monthly (22 day) volatility forecasts generated by the alternative models; for instance, the weekly volatility forecast at date T is the forecast for RK T +1 + · · · + RK T +5 (the 5-day-ahead cumulative realized kernel). Table 2 reports the percentage shares of (1, 5, 10, and 22-day) cumulative RK t out-of-sample observations that belong to the 99%, 95%, and 90% one-sided upper prediction intervals based on the distribution of the simulated sample paths; these upper prediction intervals for volatility are related to higher levels of risk in financial markets. Overall, it is seen that the empirical coverage rates of the StMAR based prediction intervals are closer to the nominal levels than the ones obtained with the reference models. By comparison, the accuracy of the prediction intervals obtained with the popular HAR model quickly degrade as the forecast period increases. The StMAR model performs well also when two-sided prediction intervals and point forecast accuracy are considered (for details, see the Supplementary Material).

Properties of the multivariate Student's t-distribution
The standard form of the density function of the multivariate Student's t-distribution with ν degrees of freedom and dimension d is (see, e.g., Kotz & Nadarajah (2004, p where Γ (·) is the gamma function and µ ∈ R d and Σ (d × d), a symmetric positive definite matrix, are parameters. For a random vector X possessing this density, the mean and covariance are E[X] = µ and Cov[X] = Γ = ν ν−2 Σ (assuming ν > 2). The density can be expressed in terms of µ and Γ as This form of the density function, denoted by t d (x; µ, Γ, ν), is used in this paper, and the notation X ∼ t d (µ, Γ, ν) is used for a random vector X possessing this density. Condition ν > 2 and positive definiteness of Γ will be tacitly assumed. For marginal and conditional distributions, partition X as X = (X 1 , X 2 ) where the components have dimensions d 1 and d 2 (d 1 + d 2 = d). Conformably partition µ and Γ as µ = (µ 1 , µ 2 ) and Then the marginal distributions of X 1 and X 2 are t d 1 (µ 1 , Γ 11 , ν) and t d 2 (µ 2 , Γ 22 , ν), respectively. The conditional distribution of X 1 given X 2 is also a t-distribution, namely (see Ding (2016, Sec. 2)) Now consider a special case: a (p + 1)-dimensional random vector X ∼ t p+1 (µ1 p+1 , Γ p+1 , ν), where µ ∈ R and Γ p+1 is a symmetric positive definite Toeplitz matrix. Note that the mean vector µ1 p+1 and the covariance matrix Γ p+1 have structures similar to those of the mean and covariance matrix of a (p + 1)-dimensional realization of a second order stationary process. More specifically, assume that Γ p+1 is the covariance matrix of a second order stationary AR(p) process.
Partition X as X = (X 1 , X 2 ) = (X 1 , X p+1 ) with X 1 and X p+1 real valued and X 1 and X 2 both p×1 vectors. The marginal distributions of X 1 and X 2 are X 1 ∼ t p (µ1 p , Γ p , ν) and X 2 ∼ t p (µ1 p , Γ p , ν), where the (symmetric positive definite Toeplitz) matrix Γ p = Cov [X 1 ] = Cov [X 2 ] is obtained from Γ p+1 by deleting the first row and first column or, equivalently, the last row and last column (here the specific structures of µ1 p+1 and Γ p+1 are used). The conditional distribution of X 1 given X 2 = x 2 is where expressions for µ(x 2 ) and σ 2 (x 2 ) can be obtained from above as follows. Partition Γ p+1 as

Supplementary material for "A mixture autoregressive model based on Student's t-distribution" by Meitz, Preve, and Saikkonen
This Supplementary Material includes proofs of Theorems 1-3, information on the numerical optimization methods employed for maximum likelihood estimation, simulation experiments, and further details of the empirical example.

Proofs
Proof of Theorem 1. Corresponding to ϕ 0 ∈ R, ϕ = (ϕ 1 , . . . , ϕ p ) ∈ S p , σ > 0, and ν > 2, define the notation Γ p , γ 0 , γ p , µ, and Γ p+1 as in (4), and note that Γ p and Γ p+1 are, by construction and due to assumption ϕ ∈ S p , symmetric positive definite Toeplitz matrices. To prove (i), we will construct a p-dimensional Markov process z t = (z t , . . . , z t−p+1 ) (t = 1, 2, . . .) with the desired properties. We need to specify an appropriate transition probability measure and an initial distribution. For the former, assume that the transition probability measure of z t is determined by the density function t 1 (z t ; µ(z t−1 ), σ 2 (z t−1 ), ν + p), where µ(z t−1 ) and σ 2 (z t−1 ) are obtained from the last two displayed equations in the Appendix by substituting z t−1 for x 2 . This shows that z t can be treated as a Markov chain (see Meyn and Tweedie (2009, Ch. 3)). Concerning the initial value z 0 , suppose it follows the t-distribution z 0 ∼ t p (µ1 p , Γ p , ν). Furthermore, if z + t = (z t , z t−1 ) = (z t , z t−p ), we find from the Appendix that the density function of z + 1 is given by Thus, z + 1 ∼ t p+1 (µ1 p+1 , Γ p+1 , ν) and, as in the Appendix, it follows that the marginal distribution of z 1 is the same as that of z 0 , that is, z 1 ∼ t p (µ1 p , Γ p , ν) (the specific structure of Γ p+1 is used here). Hence, as z t is a Markov chain, we can conclude that it has a stationary distribution characterized by the density function t p (z, µ1 p , Γ p , ν) (see Meyn and Tweedie (2009, pp. 230-231)). This completes the proof of (i).
To prove (ii), note that, due to the Markov property, z t | F z t−1 ∼ t 1 (µ(z t−1 ), σ 2 (z t−1 ), ν + p) where F z t−1 signifies the sigma-algebra generated by {z s , s < t}. Thus we can write the conditional expectation and conditional variance of z t given F z t−1 as Denote this conditional variance by σ 2 t = σ 2 (z t−1 ) (and note that σ 2 t > 0 a.s. due to the assumed conditions σ 2 > 0, Γ p > 0, and ν > 2). Now the random variables ε t defined by ε t def = (z t − ϕ 0 − ϕ ′ z t−1 )/σ t follow, conditional on F z t−1 , the t 1 (0, 1, ν + p) distribution. Hence, we obtain the 'AR(p)-ARCH(p)' representation (7). Because the conditional distribution ε t | F z t−1 ∼ t 1 (0, 1, ν + p) does not depend on F z t−1 (or, more specifically, on the random variables {z s , s < t}), the same holds true also unconditionally, ε t ∼ t 1 (0, 1, ν + p), implying that the random variables ε t are independent of F z t−1 (or of {z s , s < t}). Moreover, from the definition of the ε t 's it follows that {ε s , s < t} is a function of {z s , s < t}, and hence ε t is also independent of {ε s , s < t}. Consequently, the random variables ε t are IID t 1 (0, 1, ν + p), completing the proof of (ii).
Proof of Theorem 2. First note that y t is a Markov chain on R p . Now, let y 0 = (y 0 , . . . , y −p+1 ) be a random vector whose distribution has the density f (y 0 ; θ) = M m=1 α m t p (y 0 ; µ m 1 p , Γ m,p , ν m ). According to (8), (9), (11), and (A1), the conditional density of y 1 given y 0 is It follows that the density of (y 1 , y 0 ) is f ((y 1 , y 0 ); θ) = M m=1 α m t p+1 ((y 1 , y 0 ); µ m 1 p+1 , Γ m,p+1 , ν m ). Integrating y −p+1 out (and using the properties of marginal distributions of a multivariate t-distribution in the Appendix) shows that the density of y 1 is f (y 1 ; θ) = M m=1 α m t p (y 1 ; µ m 1 p , Γ m,p , ν m ). Therefore, y 0 and y 1 are identically distributed. As {y t } ∞ t=1 is a (time homogeneous) Markov chain, it follows that {y t } ∞ t=1 has a stationary distribution π y (·), say, characterized by the density f (·; θ) = M m=1 α m t p (·; µ m 1 p , Γ m,p , ν m ) (cf. Meyn and Tweedie (2009, pp. 230-231)). For ergodicity, let P p y (y, ·) = Pr(y p | y 0 = y) signify the p-step transition probability measure of y t . It is straightforward to check that P p y (y, ·) has a density given by The last expression makes clear that f (y p | y 0 ; θ) > 0 for all y p ∈ R p and all y 0 ∈ R p . Now, one can complete the proof that y t is ergodic in the sense of Meyn and Tweedie (2009, Ch. 13) by using arguments identical to those used in the proof of Theorem 1 in .
Proof of Theorem 3. First note that Assumption 1 together with the continuity of L (c) T (θ) ensures the existence of a measurable maximizerθ T . For strong consistency, it suffices to show that a certain uniform convergence condition and a certain identification condition hold. Specifically, the former required condition is that the conditional log-likelihood function obeys a uniform strong law of large numbers, that is, sup θ∈Θ |L , condition E [sup θ∈Θ |l t (θ)|] < ∞ ensures that the uniform law of large numbers in Ranga Rao (1962) applies.
Given conditions (i)-(iii) of the theorem, asymptotic normality of the ML estimator can now be established using standard arguments. The required steps can be found, for instance, in Kalliovirta et al. (2016, proof of Theorem 3). We omit the details for brevity.

Numerical optimization
Finding maximum likelihood estimates of the unknown parameters of an StMAR(p,M ) model amounts to maximizing L (c) T (θ), a function in M (p + 4) − 1 variables, under several constraints. Our experience with both actual and simulated data indicates that this can be challenging, in part due to multiple local maxima, and that advanced numerical optimization methods are needed. We use a hybrid numerical optimization scheme combining randomized search methods and classical gradient based methods to efficiently search for a global maximum that satisfies the constraints.
We first employ a genetic algorithm using a variety of initial populations (collections of starting points; for discussions on the genetic algorithm, other popular algorithms, and their applications in econometrics, see Goffe et al., 1994, andDorsey andMayer, 1995). For each of the initial populations, the genetic algorithm is run for a small number of generations to reach the region near an optimum point relatively quickly. Corresponding to each initial population, the solution found by the genetic algorithm is then used as a starting point for Matlab's optimization method fmincon, which is faster and more efficient for local search (for fmincon we further use a sequential quadratic programming method; see e.g. Nocedal and Wright, 2006). The final parameter estimate is the best solution found by fmincon for all the starting points considered. This hybrid optimization scheme combining multiple initial populations, the genetic algorithm, and fmincon allows us to efficiently search the parameter space and reduces the risk of ending up with a local, not global, maximum. We parallelize our code to consider multiple initial populations and starting points in parallel. This helps to speed up the optimization considerably. In view of the complexity of the estimation procedure, numerical gradients and Hessians are used for the optimization.
The StMAR code used in our S&P 500 realized kernel example, further described in our StMAR MATLAB Toolbox Documentation, is available for download through the second authors webpage at https://www.researchgate.net/profile/Daniel_Preve. R code by Savi Virolainen is available through the CRAN repository in the form of the 'uGMAR' package.

Simulation experiments
We carried out several Monte Carlo studies to evaluate the performance of the numerical optimization scheme described above. The results of two of these studies are reported in Tables 1 and 2. In these experiments, 500 independent simulated sample paths were generated from an StMAR(1,2), and also from an StMAR(4,2), process; the sample sizes and parameter values used are displayed in Tables 1  and 2.
Overall, the performance of the numerical optimization scheme is quite satisfactory. As is commonly known, the degrees of freedom parameter of a Student's t-distribution is relatively difficult to estimate, especially if its true value is large. This is also the case for our StMAR model, and our simulation results indicate that the ν m parameters can be relatively difficult to estimate even in moderate or large samples. Similar difficulties were reported by  when estimating their (constant mixing weights) version of a Student t-mixture autoregressive model using the EM algorithm (see their  Table 3).

In-sample results
We estimated 12 different StMAR models with p = 1, 2, 3, 4 and M = 1, 2, 3. Of these models, the bic, hqc, and aic information criteria chose the StMAR(4,1), StMAR(4,2), and StMAR(4,3) models, respectively. Estimation results for these three models are shown in Table 1 of the main paper. Higherorder models were also tried but their forecasting performance was inferior to the models with p = 4. Table 2 of the main paper reported the percentage shares of 1, 5, 10, and 22-day cumulative RK t out-of-sample observations that belong to the 99%, 95%, and 90% one-sided upper prediction intervals based on the distribution of the simulated sample paths. The corresponding numbers for two-sided prediction intervals (for nominal levels 99%, 95%, 90%, 70%, and 50%) are presented in Table 3.

Two-sided prediction intervals
Overall, it is seen that the empirical coverage rates of the StMAR based prediction intervals are closer to the nominal levels than the ones obtained with the reference models. The StMAR(4,1), and also the StMAR(4,2), does particularly well. By comparison, the accuracy of the prediction intervals obtained with the HAR quickly degrade as the forecast period increases. Note that to generate prediction intervals for the reference AR and HAR models, we need to specify an error distribution in these models; we assume that the errors are Gaussian. The order of the AR model is chosen using aic and bic; both favour an AR(p) model with p = 11.

Volatility point forecast evaluation criteria
Let RM denote a (cumulative) realized measure (volatility proxy), such as the realized variance or realized kernel, and RM a forecast of RM . Although realized measures generally are consistent estimators of the underlying latent volatility, in practice they are noisy proxies. Because of this, care needs to be taken when choosing a loss function to evaluate and compare volatility forecasts. Following the literature on volatility forecast comparison (Patton and Sheppard, 2009;Patton, 2011), we consider the two most widely used loss functions, namely squared loss ( Moreover, as Patton and Sheppard (2009) recommend the use of QLIKE rather than MSE in volatility forecasting applications, we employ QLIKE loss as our primary loss function, and squared loss as our secondary loss function.

Point forecasts
Results for 1, 5, 10, and 22-day cumulative RK t forecasts based on the sample median are presented in Figure 1. The left panel reports QLIKEs and the right panel MSEs. Forecast accuracy of the models is reported relative to the StMAR(4,2) model: The horizontal line (at 100) represents the StMAR(4,2) model, whereas the other lines represent the size of the forecast error measure made relative to this model (for instance, a value of 110 in the left panel is to be interpreted as a QLIKE 10% larger than for the StMAR(4,2) model). The overall performance of the StMAR(4,2) model is quite reasonable. The model does particularly well in terms of our primary loss function, QLIKE. Overall, the StMAR(4,3) performs somewhat more poorly in terms of MSE. Figure 1 also suggests that the more parsimonious StMAR(4,1) model may be preferred to the StMAR(4,2) model over longer (biweekly, monthly) forecast periods. The popular HAR model performs well under MSE, but considerably less so under QLIKE.  Figure 1: Relative forecast accuracies for the S&P 500 RK data in terms of QLIKE (left) and MSE (right). Results for the AR(11) (circle), HAR (square), StMAR(4,1) (diamond), StMAR(4,2) (solid), and StMAR(4,3) (triangle) models.