Title: MULTISCALE HIDDEN MARKOV MODELS FOR CO-VARIANCE PREDICTION

Abstract: This paper presents a novel variant of hierarchical hidden Markov models (HMMs), the multiscale hidden Markov model (MSHMM), and an associated spectral estimation and prediction scheme that is consistent, finds global optima, and is computationally efficient. Our MSHMM is a generative model of multiple HMMs evolving at different rates where the observation is a result of the additive emissions of the HMMs. While estimation is relatively straightforward, prediction for the MSHMM poses a unique challenge, which we address in this paper. Further, we show that spectral estimation of the MSHMM outperforms standard methods of predicting the asset covariance of stock prices, a widely addressed problem that is multiscale, non-stationary, and requires processing huge amounts of data.

1. INTRODUCTION:

Hidden Markov models (HMMs) are a critical tool in many fields, including speech recognition (Huang et al., 1990) , robotics (Thrun et al., 1998) , and natural language processing (NLP) (Rodu et al., 2013) . Spectral estimation of HMMs, used in place of the expectation maximization (EM) algorithm or Gibbs sampling, now allows for fast estimation, making HMMs feasible for large datasets (Hsu et al., 2012; Foster et al., 2012) . However, modeling complex structure with a standard HMM, if even possible, can require a prohibitively large state space. Previous work using HMMs to predict volatility fail to capture multiple timescales, and additionally are not able to scale to large datasets because they are estimated by EM (Rossi & Gallo, 2006; Nystrup et al., 2016) . To specifically address datasets with multiple timescales, we propose the multiscale hidden Markov model (MSHMM), which utilizes computationally efficient spectral estimation to accommodate large, high-dimensional datasets, and can be usefully applied to multiscale non-stationary problems of industrial scale. However, variants of FHMMs have been proposed for multiscale tasks such as speech recognition, natural language modeling, and energy optimization Cetin & Ostendorf (2004) ; Kolter & Jaakkola (2012) ; Aarno & Kragic (2007) ; Nepal & Yates (2013) ; Oliver et al. (2004) .

In order to motivate our method we apply the MSHMM to high-frequency asset pricing data. For US equities, sampling once per second yields roughly 5 million observations per year per stock, on thousands of stocks. Prediction of asset covariance has been recognized as an extremely important task (Wu & Xiao, 2012) , and continues to be of interest (Nystrup et al., 2017) , but modeling at that scale has proven challenging and often intractable. Since asset covariance is a non-stationary stochastic process, the MSHMM structure is well suited for this problem.

Our main contribution is a novel variant of multiscale hidden Markov model that uses a unique, structured hierarchical form to force the model to evolve at multiple timescales. The prediction phase of MSHMM differs markedly from other hierarchical HMMs due to the additive emission and multiscale noise structure. The resulting MSHMM provides a model that allows for spectral estimation, thus avoiding computationally costly EM or Gibbs sampling. Spectral estimation's ability to offer additional speed, while avoiding local minima, allows us to use our method to model large-scale multivariate time series such as high-frequency stock trades. We prove that MSHMM is consistent under spectral estimation, and we provide an efficient estimation algorithm. We show that MSHMM can achieve excellent performance in the real world task of predicting realized covariance and compare to existing methods.

2. MULTISCALE HIDDEN MARKOV MODEL:

A hidden Markov model characterizes a sequence of discrete-time observations (in this paper we consider high-dimensional, continuously distributed observations) as being emitted by a chain of discrete hidden states. There are two main assumptions, 1) that observations are conditionally independent given the hidden states, and 2) that the hidden states are governed by a Markov process, meaning that the distribution over hidden states at time t given the entire latent state history depends only on the previous hidden state at time t -1, so

Pr[h t | h t-1 , . . . , h 1 ] = Pr[h t | h t-1 ].

We can fully specify the HMM by a transition matrix T ∈ R m×m where m is the number of hidden states, an emission probability λ that specifies the conditional probability of an observation x given a hidden state h, and an initial distribution π over hidden states.

π i = Pr[h 1 = i], T i,j = Pr[h t+1 = i | h t = j], [λ(x)] i = Pr[x t |h t = i].

It is convenient to think of hidden states as being an indicator vector in R m , and λ as associating with each observation x a vector λ(x) where λ(x

) i = Pr[x|h = i].

For readers familiar with HMMs, x is often discrete, allowing λ to be easily interpreted as a matrix with finitely many rows. As previously mentioned, we assume our observations are high-dimensional x ∈ R n , and we assume that m n.

HMMs can model many behaviors exhibited by data. Flexibility is obtained by increasing the size of the hidden state space, at the cost of an often prohibitive increase in sample complexity. This can be avoided by factoring the hidden state space into multiple layers and placing a nontrivial dependence structure on these layers and emitted observations (Ghahramani & Jordan, 1997) , known as a factorial HMM (FHMM). Figure 1 (a) displays the dependence structure of the original FHMM. (2) t

h (2) t-1 h (1) t-1 Xt-1 h (2) t h (1) t Xt h (2) t+1 h (1) t+1 Xt+1 h (2) t-1 h (1) t-1 X (1) t-1 h (1) t X (1) t h (2) t+1 h (1) t+1 X (1) t+1 X (2) t-1 X (2) t+1 + + + Xt-1 Xt Xt+1 ∆ (2) = 2

is the hidden state sequence of the slower moving HMM with transitions happening at time steps ∆ (2) = 2. The observations X t are the sum of the emissions of the two HMMs, so

X t = X (1) t + X (2) t-1 .

In contrast to the FHMM, MSHMM has multiple, independently evolving HMMs with additive emissions as seen in Figure 1 (b). Consider a system with M hidden state chains {h (1) , . . . , h (M ) }. Each hidden state chain h (i) transitions at different time periods ∆ (i) , so that at a given time t, chains transition by time t + 1 only if they satisfy (t + 1 modulo ∆ (i) ) = 0, with self-transitions allowed. In this paper, we assume that the parameters {∆ (1) , . . . , ∆ (M ) } are unique, and thus without loss of generality ordered from smallest to largest, so that the hidden state chain h (M ) is the slowest evolving chain. Unlike in Ghahramani & Jordan (1997) , where observations are sampled from a Gaussian centered at a linear combination of the hidden states of each chain, in our model each hidden state chain emits an unobserved output with noise, and the system-wide observations are linear combinations of the output. In other words, we can think of each layer as an actual HMM with unobserved, though otherwise typical, output. Critically, we assume the noise terms are resampled only when there is a (possibly self-) transition.

2.1. SPECIFICATION OF MULTISCALE HIDDEN MARKOV MODEL:

Recalling that an HMM can be fully specified with parameters π, T , and λ(x), we denote the HMM parameters of each hidden process i with a superscript (i), so π (i) , T (i) , and λ (i) (x). As above, h (i) represents the hidden state chain for hidden process i, and h (i)

t its value at time t. Furthermore, m (i) denotes the size of the hidden state space for process i. Letting h t without the superscript i denote the global hidden state of the entire system, we can characterize the independence between hidden states as

Pr[h t+1 | h t ] = M i=1 Pr[h (i) t+1 | h (i) t ].

Denote the expectation of the continuous emission for model i at time t, q(h

(i) t ) = E[X (i) t |h (i) t ]

and let q(h t ) be a matrix collecting q(h

(i) t ) on the columns, so q(h t ) = q(h (1) t ), q(h (2) t ), • • • , q(h (M ) t ) .

To incorporate noise at each level into our model, we define a matrix q(h t ) as

q(h t ) = q(h (1) t ) + (1) t , • • • , q(h (M ) t ) + (M ) t

, where (i) t ∼ N (0, σ (i) ); however, the model is robust to other distributions. The observable emission is then X t = q(h t )w where w is a vector of size M . Figure 1 displays the representation of MSHMM. Note that, as a consequence of the stickiness property in our model,

Pr[q(h (i) t+1 ) + (i) t+1 = q(h (i) t ) + (i) t ] = 1 if (t + 1 modulo ∆ (i) ) = 0. Finally, we let h(i) t = E[h (i) t |X t-1 . . . X 1 ]

denote the belief state of HMM i at time t.

2.2. PREDICTION MULTISCALE HIDDEN MARKOV MODEL:

The prediction procedure for our model differs substantially from standard hierarchical and factorial HMMs. Typically for HMMs and variants of HMMs, it is only necessary to estimate the belief state ht . But for the MSHMM, the correction term

˜ (i) t q(h (i) t ) + (i)

t -q( h(i) t ) from HMM i acts as persistent bias for HMM i-1 for a period ∆ (i) . Therefore, in order to eliminate this bias, we must estimate

˜ (i) t .

We require one more definition in order to specify the prediction procedure. Let g (i) (t) to be a step function for the ith hidden state chain that steps when HMM process i has a transition. For example, if t = {0, 1, 2, 3, 4, 5, 6, • • • }, and ∆ (3) = 3 then g (3) (t) = {0, 0, 0, 3, 3, 3, 6, • • • } which corresponds to when process 3 has a state transition.

Prediction proceeds in two steps. First, we obtain predictions q( h(i) t ) from each component HMM. Second, in order to eliminate the persistent bias effect mentioned above, we estimate the noise terms ˜ (i) t for i = 2 . . . M and add these to q( h(i) t ) to obtain an estimate of q(h

(i) t ) + (i) t . We first estimate ˜ (M )

t+1 , where M is the slowest HMM process, as

E[˜ (M ) t |h t ,X t , • • • , X g (M ) (t) ] = 1 t -g (M ) (t) t-g (M ) (t) j=0 X t-j -q( h(M) t )w (M ) ,

and subsequently estimate

E[˜ (i) t |h t ,X t , • • • , X t-g (i) (t) ] = 1 t -g (i) (t) t-g (i) (t) j=0 X (i) resid,t-j -q( h(i) t )w (i) ,

and finally predict

X t+1 = q(h t )w = q( h(1) t ) + ˜ (1) t , • • • , q( h(M) t ) + ˜ (M ) t w.

This correction term estimation procedure must be repeated at all times t for all component HMMs. Algorithm 1 details the prediction procedure.

Algorithm 1 Prediction using MSHMM

1: Input: Observed emissions X. 2: Let X M resid = X, and set i = M . 3: while i > 0 do 4: Predict X (i)

t+1 using spectral estimation HMM, q( h(i) t ) (see Rodu et al. (2013) for details) from the observations X (M ) resid .

5::

Compute the error

X -q( h(i) t )w (i) over (g M (t), t]. 6: Estimate new ˜ (i)

t by averaging the error over this period. 7:

Compute

X (i-1) resid = X (i) resid -q( h(i) t )w (i) -˜ (i)

t and set i = i -1. 8: end while 9: Sum over all predictions from each process X (i) t+1 .

2.3. ESTIMATION ALGORITHM FOR MULTISCALE HIDDEN MARKOV MODELS:

We divide the algorithm into two components. The first estimates the HMM parameters through the spectral method of moments (see Rodu (2014) for details). The second, presented in algorithm 2, estimates the weight vector that combines the emissions of the component HMMs in the MSHMM.

The critical assumptions that allow our model to be estimated using fast and consistent spectral methods are 1) the independence of the component HMM processes and 2) the fact that the emissions X t are linear combinations of the latent emissions from the component HMM processes. Note that the decoupling of the HMMs is possible because they evolve at different rates.

Algorithm 2 Computing weights for the linear combination of HMMs 1: Input: Training set of observed emissions X. 2: Let X (M ) resid = X, and set i = M . 3: while i > 0 do 4:

Using spectral estimation method from Foster et al. ( 2012) (see supplementary materials), estimate the HMM process i with a binned average

X (i) avg(∆ (i) )

of bin width ∆ (i) over time series

X (i)

resid , and compute observation expectations q(h (i) t ).

In order to recenter, compute

X (i-1) resid = X (i) resid -X (i) avg(∆ (i) )

and set i = i -1. 6: end while 7: Estimate w with equation Xt = q(ht)w + t with q(ht) = q(h

(1) t ), q(h (2) t ), . . . , q(h (M ) t

) .

2.4. THEORETICAL FOUNDATIONS OF MULTISCALE HIDDEN MARKOV MODELS:

In our estimation scheme, observations for HMM process M are obtained by segmenting the time series into bins of size ∆ (M ) . A residual time series is then constructed by subtracting the derived observations from their respective observations. We then obtain observations for HMM process M -1 by binning the residual time series in a similar fashion, and proceed recursively for every other process. Each time we bin, however, there is noise in the estimated observation, and that noise gets propagated to faster chains. Critically, the only effect of this noise on consistency is bias in the estimation of the diagonal elements of the second moment, which only appears in the method of moments estimation in its inverted form.

More technically, note that in our estimation scheme we assume that E[X (i) t ] = 0 for i < M . In an MSHMM with only two processes, for instance, this implies that E[X (1) t ] = 0. Let q (1) be a matrix such that

q (1) h (1) t = q(h (1) t ) = E[X (1) t |h (1) t ]

then under stationarity conditions we have that q (1) π (1) = 0. This implies that E

[X t ] = X (2) t for t ∈ [t, t + ∆ (2)

), which suggests that we first estimate observations X

(2) t+k∆ (2) of the slow-moving HMM by calculating the expected value of X t for t ∈ [t, t + ∆ (2) ), then estimating the slowmoving model using these calculated observations and the fast-moving model using the residuals of observations from their expected value. But because we have only S = ∆ (2) /∆ (1) observations X t for t ∈ [t, t+∆ (2) ) our estimate is noisy. The main problem is that the bias in the empirical estimate of

X (2) t is reflected in X t -X (2) t for all t ∈ [t, t + ∆ (2)

). This, in turn, biases the estimation of the diagonal elements of the second moment matrix of the faster HMM 1. It does not, however, bias estimation of HMM 2, as the error between two successive observations for that HMM are independent.

In practice, the bias seems to have no real effect on the performance of the MSHMM, and even might serve as a natural regularization. It is possible, however, to estimate the spectral parameters so as to eliminate this bias, which we lay out in theorem 1. This provides a practical mechanism for obtaining consistent estimates that effectively restricts the sample size of each component HMM process to be on the order of that of the slowest HMM. On the other hand, if we let the ratio between sampling intervals become arbitrarily large as in theorem 2, while impractical, this allows observations estimates from each HMM component to become arbitrarily precise, thus again providing consistent estimators. Proof. We consider the two layer process, though this generalizes to an arbitrary number of layers. As a reminder, we can decompose the theoretical contributions from the faster-moving HMM 1 into two components: the state-conditional means and a mean 0 noise term that we consider to be drawn iid at each time step. Similarly, empirical contributions can be decomposed as the state-conditional means, an iid noise term, and an estimation bias propagated from HMM 2 that persist for ∆ (2) time steps. Note that from the central limit theorem the bias terms are drawn iid with mean 0. Therefore when estimating the second order observable as suggested in theorem 1 the bias term simply acts as an additional noise term. It is straightforward to see, then, that all the necessary observables converge to their true parameters.

Theorem 2. Assume we have an MSHMM as defined above, and let J be the total number of observations. Let J → ∞ and ∆ (i) /∆ (i-1) → ∞ for 1 < i ≤ M . Then the estimated spectral MSHMM parameters converge to the true spectral MSHMM parameters (in the sense that the probability forecasts are consistently estimated.)

Proof. Considering the two layer MSHMM, we show that as

EQUATION

t → X

(2)

t . We write

X(2) t = 1 S t ∈[t,t+∆ (2) ) X t = X (2) t + q (1) π + 1 S t ∈(t,t+∆ (2) )

(1) t

From the law of large numbers, as S → ∞, 1 S t ∈(t,t+∆ (2) )

(1) t → 0, and π(1) → π (1) , and since

q (1) π (1) = 0 we have X(2) t → X (2) t .

While theorem 2 is impractical, it shows that as the separation of time scales between HMM processes increases, the empirical estimates of the spectral parameters converge to their theoretical counterparts. If a comfortably large enough separation is not feasible or desirable, theorem 1 provides an estimation method that inherently enjoys the usual consistency of spectral HMM estimation.

In practice, however, the effect of the bias is minimal.

3 EXPERIMENTS WITH SYNTHETIC DATA We assessed the estimation speed and accuracy of MSHMMs using synthetic data generated from an MSHMM process containing 3 component HMMs each with 5 hidden states (MSHMM-3-5) which emits a vector of continuous observations. Computation time for training the MSHMM was linear in the number of observations and quadratic in the number of states. For our application, we found no benefit beyond 15 hidden states. Our computational times are consistent with Parikh et al. ( 2012). Furthermore, because of the MSHMM specification, the HMM with the highest resolution dominated the computational time for estimating the model. For comparison, a simple HMM with 5 hidden states using EM1 required 1255 seconds to estimate parameters for 900,000 observations while our MSHMM-3-5 took 25 seconds.

In Figure 2 we show our results using synthetic data. We compared the correct MSHMM to HMMs with 5 and 15 hidden states. We also compared a 3 Markov component MSHMM which has 10 hidden states instead of 5. Finally, as an upper bound on performance we used the "oracle" MSHMM which makes prediction using the true hidden state transition, emission, and weight parameters, where only the hidden state sequence is unknown. Surprisingly, the MSHMM with too many hidden states still performs well. The relative root-mean-squared error (RMSE) is the RMSE normalized by the oracle RMSE. This points to the fact that the model is robust to some misspecification; however, further regularization might improve results as shown by Tran et al. (2016) .

We also compared our model against Long-Short Term Memory (LSTM) (Hochreiter & Schmidhuber, 1997) recurrent networks, which with a 30 dimensional hidden state trained using all 200,000 data points gave a relative RMSE of 1.76. While this is much better than the results from HMM-15, the LSTM failed to capture multiscale long term dependencies and so underperformed our MSHMM. We further examined the State Frequency Memory (SFM) recurrent network (Hu & Qi, 2017) , which incorporates frequency information into the memory cells of an LSTM. Its relative RSME was 1.69 using 10 frequency banks, again giving poorer accuracy than the MSHMM

4. ESTIMATING STOCK PRICE COVARIANCE:

4.1 DATA DESCRIPTION Daily stock closing prices were extracted from the Center for Research in Securities Prices (CRSP) database.2 All intraday data was extracted from the NYSE Trade and Quote (TAQ) database. 3 All stock returns are calculated using the difference of log prices. The calculation of realized covariance uses the same procedure as Sheppard & Xu (2014) .

We sampled the realized covariance at daily, 17-minute and 1-minute rates. We chose the 17-minute sample rate since Bollerslev et al. (2008) have shown that microstructure effects on realized variance are mostly mitigated. We chose six exchange-traded funds (ETFs), including SPY (which tracks the S&P 500 index), TLT, FEZ, USO, UNG, and GLD. Each realized covariance matrix has 21 unique entries. We calculated realized variance over 10 years of historical daily data, 2 years of 17-minute intraday data, and 1 year of 1-minute intraday data. To be exact, for training we used 2268 daily, ∼ 22000 17-minute, and ∼ 20000 1-minute observations, and for testing we used 252 daily, ∼ 2600 17-minute, and ∼ 4000 1-minute observations. Due to market conditions such as stock halts, the number of intraday data points is not the same for all stocks, so we reported approximate number of observations.

4.2. EMPIRICAL RESULTS:

Our model was applied to high-frequency equity data for covariance estimation. We estimated the slowest process using 9 years of data, while the minute-level process was estimated using shorter periods of time. This is another useful feature of our model versus other hierarchical HMMs: we can estimate the daily covariance process using a longer overlapping time period.

We used the standard Generalized Autoregressive Conditional Heteroskedastic model (GARCH) (Bollerslev, 1986) as our baseline. We chose the baseline model to be GARCH(1,1), such that

σ t = h t v t , v t ∼ N (0, 1) and h t = ζ + ασ 2 t-1 + γh t-1 .

The prediction from GARCH(1,1) is the estimate of volatility σ 2 t . GARCH is a standard baseline for volatility prediction. We also compared our model to two types of models. The first is a generalization of the multivariate GARCH, Bayesian Multivariate Dynamic Covariance model (BMDC)4 (Wu et al., 2013) . The second type of model is a generalization of regime-switching models to high-frequency data, High-Frequency-Based Volatility (HEAVY)5 models (Noureldin et al., 2012; Sheppard & Xu, 2014) . BDMC and HEAVY have been shown to be computationally efficient to estimate as well as to perform at state-of-the-art levels. The model seeks to predict the covariance either a day, 17 minutes, or a minute forward. Table 1 : The RMSE of the out of sample predicted versus realized stock covariance which is normalized by the GARCH(1,1) RMSE. Run times are in parenthesis in seconds. The training time for neural models are between 20 minutes and 5 hours using a single Titan X GPU and highly variable.

Horizon:

We report the RMSE of the models normalized with respect to the GARCH baseline. Table 1 shows the results of the prediction on ETFs. MSHMM is clearly superior at higher frequency, where the realized covariance matrix tends to be sparse. This result is consistent with our simulated data. When we used a daily prediction horizon, MSHMM underperformed BDMC, HEAVY, and LSTM models; this should be expected, however, since high-frequency data does not aid estimation very much when estimating the covariance only at the market close. MSHMM outperformed both models when predicting at 17-minute intervals and 1-minute realized covariance. MSHMM is also significantly faster than BDMC and HEAVY models6 . MSHMM outperforms the LSTM and SFM neural models7 in part due to the fact that there are daily jumps between the market close and subsequent open on the next trading day. We needed to remove the overnight return effect in order to improve upon GARCH(1,1) performance. Finally, we again draw the reader's attention to the fact that without the estimation of noise as done in algorithm 1, MSHMM underperformed even GARCH(1,1).

5. RELATED WORK:

The MSHMM is a specialized hierarchical HMM, in a similar vein proposed by El Hihi & Bengio (1995) , and is distinct from the factorial hidden Markov models (FHMMs) of Ghahramani & Jordan (1997) . While FHMMs and MSHMMs are both capable of learning short-and long-term structure, the MSHMM explicitly differentiates time-scale evolution.

Spectral estimation of latent state models Recently there has been significant interest in spectral estimation of latent state models, beginning with the seminal work in spectral estimation of discrete HMMs (Hsu et al., 2012) which has been extended to fully reduced methods (Foster et al., 2012) , to continuous output HMMs (Song et al., 2010; Rodu et al., 2013) , to trees, (Dhillon et al., 2012) , and to more general latent variable structures (Anandkumar et al., 2014; Parikh et al., 2012) . In this paper, we extend spectral estimation techniques to MSHMMs.

Hierarchical HMM variants A multitude of hierarchical HMM variants have been proposed since hierarchical HMMs were first introduced (El Hihi & Bengio, 1995) . Hierarchical HMMs are characterized by multiple interconnected hidden states, Pr[h t+1 | h t , h t , • • • ] where h t is not in the same hidden state sequence, or by multiple hidden states that affect an emission, Pr[

x t | h t , h t , • • • ]. FH-

MMs are a subclass of hierarchical HMMs in which multiple independent hidden state processes affect an observable emission. There are a multitude of variants on FHMMs. Among these are multirate HMMs (Cetin & Ostendorf, 2004 ) and additive FHMMs (Kolter & Jaakkola, 2012) , which differ in structure from our model.

Multirate HMMs have multiple hidden state processes, which we refer to as "levels" or "components," each evolving at different rates. The main structural difference between the multirate HMM and our model is that only the first level emits an observable; all subsequent "higher" levels output to the level "below," which has a faster hidden state process. MSHMM has additive emissions results, whereas, multirate HMM will propagate the same error upstream through different levels of hidden states.

By way of contrast, in additive FHMMs all hidden state processes evolve at the same rate, but the higher level process emits to both observations at t and t + 1. It is possible to formulate an MSHMM and more generally hierarchical HMM as a latent junction tree, on which spectral algorithms can be estimated (Parikh et al., 2012) . However, latent junction trees are inefficient representations for estimating FHMMs and MSHMMs (Jordan et al., 1999) as they couple the hidden state paths, leading to an explosion in the number of parameters to be estimated.

Neural Models Deep neural networks allow a non-linear form of hierarchical HMMs as seen in the work of El Hihi & Bengio (1995) and more recently Liu et al. (2015) ; Chung et al. (2016) . Similar to junction trees, these models require exponentially large hidden state representations (Jaderberg et al., 2016) compared to MSHMM. The work of Chung et al. (2016) is more similar to multi-rate HMMs than to MSHMM. Nonetheless, exploration of a neural equivalent to MSHMM would be interesting. The SFM recurrent networks (Hu & Qi, 2017; Zhang et al., 2017) have been shown to be useful for stock price prediction, which is a related task. However, SFM does not explicitly incorporate multiple time scales, which is critical for covariance estimation. Another recent neural model is the neural volitility model (Luo et al., 2017) which has shown superior results than GARCH models; however, this model is not multiscale.

Financial Covariance Estimation Models:

The standard methods for modeling stock price variation trade off between parameter accuracy and computational efficiency. The modern classes of models are multivariate versions of generalized autoregressive conditional heteroscedasticity (GARCH) (Engle & Kroner, 1995; Wu et al., 2013) , stochastic volatility models (Gouriéroux et al., 2009) , and Markov-switching processes (Calvet & Fisher, 2008) , (which, while a hidden state model, their formulation and estimation are significantly different). Furrer & Bengtsson (2007) used a Kalman filter approach in order to estimate covariance matrices. More recently, Nystrup et al. (2016) modelled covariance using HMMs. Most of these models do not scale to large amounts of data and are used to predict daily returns or even longer periods. For a more detailed overview of models for large covariance matrix estimation models see Fan et al. (2016) .

With the emergence of high-frequency data containing detailed records of quotes and transaction prices with nanosecond time resolution, a new class of models has emerged, such as multivariate high-frequency-based-models (HEAVY) (Noureldin et al., 2012) . Our model extends the work of Rossi & Gallo (2006) to multiple timescales, covariance estimation, and the use of computationally efficient spectral methods instead of EM.

6. CONCLUSION:

Our MSHMM is a specialized hierarchical HMM that captures multiscale structure by adding emissions from multiple component HMMs. Unlike standard factorial HMMs, MSHMM handles noise at multiple time scales in a manner that cleanly supports spectral estimation. This estimation is straightforward if one accounts for the noise that propagates from the slower to the faster time scales. However, prediction requires accounting for the fact that estimation error due to noise at slower time scales acts as a bias term for the faster component HMMs. We have provided fast estimation schemes with provable asymptotic properties.

Because we cast MSHMMs in the spectral estimation framework, unlike other hierarchical HMMs, they scale gracefully with the size and dimensionality of the data, increasing the range of problems they can address. They apply naturally to high-frequency financial time series, as they capture both the short-and long-term covariances between stock prices, which are used as the basis for decisionmaking about portfolio management and risk assessment, and they rapidly detect market changes.

The MSHMM can also potentially be applied to other types of multi-time scale data in domains such as natural language processing, vision, speech recognition, neuroscience, and macroeconomics.

Supplemental Material:

A SPECTRAL ESTIMATION ALGORITHM 

: Compute μ = 1 N T N T i=1 y1 i , Σ = 1 N T N T i=1 y2 i y 1 i and tensor Ĉ = 1 N T N T i=1 y3 i ⊗ y1 i ⊗ y2 i . 6: Set b1 = μ and b ∞ = b 1 Σ-1 .

7: Right multiply each slice of the tensor in the y2 direction (so y2 is being sliced up, leaving the y3y 1 matrices intact) by Σ-1 to form B(γ) = Ĉ(γ) Σ-1 .

Spectral estimation can be understood through the lens of the observable operator model (Jaeger, 2000) . For simplicity of presentation we describe spectral estimation for a discrete emission HMM. Consider an HMM where T is an m×m transition matrix on the hidden state, O is a v ×m emission matrix giving the probabilities of hidden state h = j emitting observation x = i, and π is a vector of initial state probabilities in which πi is the probability that h1 = i. Jaeger Jaeger (2000) showed that the joint probability of a sequence of observations from this HMM is given by

P r(x1, x2, . . . , xt) = 1 Ax t Ax t-1 • • • Ax 1 π,

where Ax ≡ T diag(O δx), δx is the unit vector of length v with a single 1 in the xth position and diag(v) creates a matrix with the elements of the vector v on its diagonal and zeros everywhere else. Ax is called an 'observation operator'. Now, define a random variable yt = U δx t , where U has orthonormal columns and is a matrix mapping from observations to the reduced dimension space, then and µ = E(y1), Σ = E(y2y 1 ), and C(y) = E(y3y 1 y 2 )y. Writing µ, Σ, and C(y) in terms of their theoretical quantities, we get

µ = U Oπ Σ = U O T diag(π) O U Σ -1 = (O U ) -1 diag(π) -1 T -1 (U O) -1 C(y) = U O T diag(O U y) T diag(π) O U

and thus can be shown that b1

= U O π b ∞ = 1 (U O) -1 By = U O Ax (U O) -1

and thus

P r(x1, x2, . . . , xt) = 1 Ax t Ax t-1 • • • Ax 1 π = b ∞ By t By t-1 • • • By 1 b1

B COMPARISON TO RELATED HIDDEN MARKOV MODEL Multirate HMMs have multiple hidden state processes, which we refer to as "levels" or "components," each evolving at different rates. The main structural difference between the multirate HMM (as seen in figure 1 ) and our model (figure 3 ) is that only the first level emits an observable; all subsequent "higher" levels output to the level "below," which has a faster hidden state process. MSHMM has additive emissions results, whereas, multirate HMM will propagate the same error upstream through different levels of hidden states.

EQUATION

Xt+1 h

(1) t+2

Xt+2:

Figure 3 : Multi-rate HMM with 2 hidden state processes.

By way of contrast, in additive FHMMs (figure 2 ) all hidden state processes evolve at the same rate, but the higher level process emits to both observations at t and t + 1.

EQUATION

Xt+1 

D FURTHER DISCUSSION OF OTHER HIERARCHICAL HMMS:

Further, latent junction trees are constrained to a fixed tree topology. This has two ramifications. First, a separate model needs to be estimated for each desired topology, and more importantly, in the application presented here, Daily returns where adjusted for corporate actions such as dividends and stock splits. All intraday data was extracted from the NYSE Trade and Quote (TAQ) database10 . Intraday prices where filtered using trade flags and incorrect price quotations. Our reference intraday prices used the average of the best bid price and best ask price, and all stock returns are calculated using the difference of log prices.

Figure 6 is illustrative of the multiple scales of stock variance of SPY. In particular, note that there are clearly different states for daily and intraday realized variance.

We chose the 30 stocks: TROW, AAPL, AA, CMA, FITB, BEN, MMC, HRB, BK, AFL, AXP, CB, IBM, MSFT, PNC, XOM, PGR, KEY, LM, SLM, AIG, STI, ALL, COF, UNH, SPY, GLD, TLT, FEZ, USO, UNG , as well as exchange trader funds (ETFs) such as SPY, which tracks the S&P 500 index, and GLD which tracks the value of gold prices. Each realized covariance matrix has 465 unique entries. We calculated realized variance over 10 years of historical daily data, 2 years of 17 minute intraday data, and 1 year of minute intraday data. To be exact, for training we used 2268 daily, ∼ 22000 17 minute, and ∼ 20000 minute datapoints, and for testing we used 252 daily, ∼ 2600 17 minute, and ∼ 4000 minute datapoints. Due to market conditions such as stock halts, the number of intraday data points is not the same for all stocks, so we reported approximate number of observations. Particularly, the top subplot of figure 7 corresponds to the adjusted daily prices of SPY and GLD.

Figure 7 also demonstrates the covariance between SPY and GLD, illustrating the multiscale nature of the covariance between SPY and GLD.

F FURTHER DISCUSSION OF EVALUATION

In the case of covariance prediction we define

M SE = 1 T -k T -k 1 (C -E(C))) 2 .

A close analysis of the results shows that HEAVY performs better at estimating covariance of stocks but not ETFs (SPY, GLD, USO, etc ). We surmised that this is because the covariance matrix for stocks tends to be close to a diagonal matrix, unlike ETFs, which have more significant off diagonal components. Since the variance is larger than the covariance, this accounted for the outperformance. We surmised that this is because companies are driven by very few factors, as a result the MSHMM hidden states focused on predicting covariance of ETF. However, at relatively high frequency timefames, our model outperformed. Again the intuition is that microstructure creates both fat tails as well as sparse covariance matrices, which are akin to discrete problems. 

:

http://cran.r-project.org/web/packages/HMMCont/HMMCont.pdf

http://www.crsp.com/products/research-products/crsp-us-stock-databases

http://www.nyxdata.com/doc/2549

Code is from https://bitbucket.org/jmh233/bmdc_icml2013

Code is from https://www.kevinsheppard.com/MFE_Toolbox

We used the highly optimized C++ code of HEAVY and the Matlab code from BDMC. Our model is written in R and not yet optimized for performance.

Code is from https://github.com/z331565360/State-Frequency-Memory-stock-prediction

http://www.crsp.com/products/research-products/crsp-us-stock-databases

http://www.nyxdata.com/doc/2549

Figures and Tables:

Figure 1: (a) Factorial HMM with 2 hidden state layers. (b) An MSHMM with two HMMs where h

Theorem 1. Assume we have an MSHMM. Let k ∆ (M ) i and k ∆ (M ) i + 1 be the last observation in HMM k before the i th state transition of HMM M and the first observation in HMM k immediately following the i th transition of HMM M , respectively. Note that for HMM k the emission value at time t = ∆ (M ) i is equal to the previous ∆ (k) time steps. Assume that the second order observables for all chains k ∈ {1, . . . , M } are estimated only using the bigram k ∆ (M ) i and k ∆ (M ) i + 1. Then the estimated spectral MSHMM parameters converge to the true spectral MSHMM parameters.

Figure 2: Model RMSE for MSHMM normalized by the RMSE of the oracle. Thus 1 is the absolute minimum.

Computing observables for spectral estimation of an HMM, fully reduced third moment1: Input: Training examples-{x1 i , x2 i , x3 i } for i ∈ {1, . . . , N T }. 2: Compute Ê[x2 ⊗ x1] x2 i x 1 i . 3:Compute 8 the left k singular vectors corresponding to the top k singular values of Ê[x2 ⊗ x1]. Call the matrix of these singular vectors Û . 4: Reduce data: ŷ = Û x. 5

P r(x1, x2, . . . , xt) = b ∞ By t By t-1 • • • By 1 b1 holds where b1 = µ b ∞ = µ Σ -1 By = C(y)Σ -1

Figure 4: Factorial HMM with M hidden state processes.

Figure 7: Computational time for MSHMM

Figure 8: S&P 500 ETF (SPY) realized variance at daily, 17-minute, and 1-minute sampling rates over 10 years, 1 year, and 1 month respectively.

Figure 9: S&P 500 ETF (SPY) and SPDR Gold Shares ETF (GLD) daily prices, and their realized daily covariance.

Peer review (iclr2018):
 The paper focuses on a very particular HMM structure which involves multiple, independent HMMs. Each HMM emits an unobserved output with an explicit duration period. This explicit duration modelling captures multiple scale of temporal resolution. The actual observations are a weighted linear combination of the emissions from each latent HMM. The structure allows for fast inference using a spectral approach. I found the paper unclear and lacking in detail in several key aspects: 1. It is unclear to me from Algorithm 2 how the weight vectors w are estimated. This is not adequately explained in the section on estimation. 2. The authors make the assumption that each HMM injects noise into the unobserved output which then gets propagated into the overall observation. What are reasons for his choice of model over a simpler model where the output of each HMM is uncorrupted? 3. The simulation example does not really demonstrate the ability of the MSHMM to do anything other than recover structure from data simulated under an MSHMM. It would be more interesting to apply to data simulated under non-Markovian or other setups that would enable richer frequency structures to be included and the ability of MSHMM to capture these. 4. The real data experiments shows some improvements in predictive accuracy with fast inference. However, the authors do not give a sufficiently broad exploration of the representations learnt by the model which allows us to understand the regimes in which the model would be advantageous. Overall, the paper presents an interesting approach but the work lacks maturity. Furthermore, simulation and real data examples to explore the properties and utility of the method are required.