*Article* **Bivariate Volatility Modeling with High-Frequency Data**

#### **Marius Matei 1,2,3,\*, Xari Rovira <sup>4</sup> and Núria Agell <sup>4</sup>**


Received: 6 August 2018; Accepted: 6 September 2019; Published: 15 September 2019

**Abstract:** We propose a methodology to include night volatility estimates in the day volatility modeling problem with high-frequency data in a realized generalized autoregressive conditional heteroskedasticity (GARCH) framework, which takes advantage of the natural relationship between the realized measure and the conditional variance. This improves volatility modeling by adding, in a two-factor structure, information on latent processes that occur while markets are closed but captures the leverage effect and maintains a mathematical structure that facilitates volatility estimation. A class of bivariate models that includes intraday, day, and night volatility estimates is proposed and was empirically tested to confirm whether using night volatility information improves the day volatility estimation. The results indicate a forecasting improvement using bivariate models over those that do not include night volatility estimates.

**Keywords:** high-frequency; volatility; forecasting; realized measures; bivariate GARCH

**JEL Classification:** C32; C53; C58

#### **1. Introduction**

We aim to improve volatility modeling by adding information that exists on latent volatility processes while the markets are closed and no transactions occur. We build upon the observation that the price at market closing usually differs from the price at market opening, despite no transactions occurring between the two recordings. Models previously proposed usually estimate volatility by including information on past day and intraday volatility, estimated from day-recorded prices and sampled at various time intervals. Some papers have proposed methods to address overnight returns. The latent volatility component apparent in periods when markets are closed, highlighted by the difference between the two prices, may be the effect of events that occurred during the market closing, both domestic or international, or may be due to other latent factors that usually influence the financial markets, and may prove useful in volatility modeling. We propose an estimation of this night latent volatility and suggest a new model that uses day, intraday, and night volatility information to model day volatility. What distinguishes our contribution from other papers published on similar topics is that we propose a two-factor structure in a realized generalized autoregressive conditional heteroskedasticity (GARCH) setting that takes advantage of the natural relationship between the realized measure and the conditional (day and night) variance. The mathematical structure is thus elegant, facilitates volatility estimation, and allows the inclusion of return-volatility dependence. We call the structure bivariate

because it uses both day and night volatility information, as opposed to the univariate ones that only use day information. To strengthen the robustness of our empirical research, we further extended this idea to a number of realized GARCH models that use day and intraday volatility information, creating an equivalent set of bivariate models that additionally use night volatility information. We obtained a class of realized GARCH models that incorporate day, night, and intraday volatility measures; they were assessed against their counterparts that did not include night volatility information using an extended set of 10 stock prices. Empirical results of the forecasting performance assessment show a degree of improvement of the newly proposed models over those that do not include night volatility measures. This finding suggests the potential of our method for volatility forecasting problems for financial assets and other assets with night latent volatility information.

Financial volatility modeling has benefited significantly from the availability of high-frequency data. The main interest in modeling using frequently sampled information and integrating it into models built to estimate day conditional variance was initiated by Andersen and Bollerslev (1998), who used realized volatility estimates extracted from intraday data (realized variance) as better estimates of conditional volatility than squared returns. They proved that by adding up squared intraday returns, the forecasted volatility would correlate closely to the future latent volatility factor.

Engle (2002) was among the first econometricians who extended the standard GARCH model to include an exogenous realized measure (the realized variance) in the conditional variance (GARCH) equation. In this model, the realized measures' variation is not explained; thus, such models (GARCH-X) are considered incomplete. Engle and Gallo (2006) proposed the multiplicative error model (MEM), which was the first attempt to contain a separate GARCH structure equation for the realized measure. A similar complete model nested in a MEM setting is the high frequency based volatility (HEAVY) model of Shephard and Sheppard (2010). Both MEM and HEAVY models are difficult to use as they work with multiple latent processes—for every realized measure used, there is a corresponding latent volatility process. The Realized GARCH model proposed by Hansen et al. (2012) combines a GARCH structure for returns with realized measures of volatility. Compared with MEM and HEAVY models, the Realized GARCH model takes advantage of the natural relationship between the realized measure and the conditional variance. Instead of introducing additional latent factors, it proposes a single measurement equation in which the realized measure is a consistent estimator of the integrated variance. Besides its elegant mathematical structure, the Realized GARCH model is easy to estimate, captures the return-volatility dependence (leverage effect), and has been empirically shown to outperform conventional GARCH. A more robust version of the Realized GARCH model was introduced by Banulescu-Radu et al. (2019), suggesting a variant that is less sensitive to outliers and minimizes the impact on volatility of days with extreme negative volatility shocks. A realized exponential GARCH model that can use multiple realized volatility measures for the modeling of a return series, using a similar framework, has also been proposed (Hansen and Huang 2016). Finding that the Realized GARCH model was insufficient for capturing the long memory of underlying volatility, Huang et al. (2016) developed a parsimonious variant of the Realized GARCH model by introducing Corsi's (2009) heterogeneous autoregressive (HAR) specification in the volatility dynamics. A multivariate GARCH model that incorporates realized measures of variances and covariances was also introduced by Hansen et al. (2014), but it did not suggest the introduction of night volatility information. Bollerslev et al. (2018) proposed asymmetric multivariate volatility models that exploit estimates of variances and covariances based on the signs of high-frequency returns to allow for more nuanced responses to positive and negative return shocks than the threshold leverage effect. Hansen et al. (2019) proposed a multivariate GARCH model that incorporates realized measures for the covariance matrix of returns.

Overnight (close-to-open) volatility is usually higher than the five-minute realized volatility estimated during trading hours, and the close-to-open price differential may trigger a distorting effect on the realized volatility. Thus, the inclusion of overnight returns when constructing the realized conditional covariance matrix of the daily returns has been empirically documented to reduce

information loss and consequently improve volatility forecasting. A common approach to account for volatility during the market's closing hours has been to calculate a close-to-open return from the price change recorded between the trading day closing and the next trading day opening, and then add its squared value to the sum of intraday returns (Bollerslev et al. 2009; Martens 2002; Blair et al. 2001). Hansen and Lunde (2005) compounded optimal weights corresponding to overnight returns and to the sum of intraday returns, and Fleming and Kirby (2011) and Fuertes and Olmo (2013) further applied it. De Pooter et al. (2008) and Fleming et al. (2003) computed it in matrix form by incorporating the cross-product of the vector of overnight returns in the summation of the matrix that provided the covariance matrix of the daily returns, acknowledging that the outer product of the vector of overnight returns is an inaccurate estimator of the integrated covariance matrix for the period when markets were closed (Fleming et al. 2003). Koopman et al. (2005); Martens (2002); and Angelidis and Degiannakis (2008) excluded the noisy overnight returns to compute an estimate of volatility during trading hours, instead of daily volatility; then, they scaled up the sum of intraday returns to cover the whole 24-h day. The literature has not yet reached a consensus on the best method of accounting for overnight returns; however, Ahoniemi and Lanne (2013) suggested that the weighted sum of the squared overnight return and the sum of intraday squared returns was the most accurate measure of realized volatility for the Standard&Poor's' S&P 500 index.

This paper suggests a method of capturing and incorporating night volatility into the day conditional volatility equation of one low-frequency as well as a number of high-frequency GARCH models. We propose a two-factor structure of the conditional variance, one for night and one for day variance, in a realized GARCH setting that takes advantage of the natural relationship between the realized measure and the conditional (day and night) variance. The mathematical structure is thus elegant, facilitates volatility estimation, and allows the inclusion of the return-volatility dependence. A general framework is formulated; based on it, a set of GARCH models is adapted such that it uses the estimation of night latent volatility to model day conditional volatility. This approach enabled us to document, in an empirical context, whether the introduction of the night volatility component, in the two-factor structure and realized GARCH setting we propose, improved the volatility modeling for each of the models discussed. The new models are called bivariate as they use both night and day volatility information and are defined to work in typical financial settings, such as volatility modeling of stock and commodity prices. We assessed the performance of the bivariate models by comparing the error functions of the forecasts of the bivariate models with those obtained when the simple versions of the models, which do not use night volatility information, were used. We call the latter models univariate models. The scope of this study was thus to analyze whether the use of night volatility information in the forms proposed improves the modeling of day volatility.

The paper proceeds as follows. Section 2 proposes the new set of bivariate realized models. Section 3 describes the data and methodology, and Section 4 summarizes the results. The paper concludes with Section 5, where final remarks are presented, and some future lines of research are proposed.

#### **2. Bivariate Realized Models**

#### *2.1. Base Model*

Existing high-frequency GARCH models estimate day conditional variance using day and intraday volatility information. We developed a class of realized models that allow constructing day volatility estimates with day, intraday, and night volatility information. Models previously proposed use return and volatility information estimated from trades that occurred during the trading day to estimate next-day volatility. However, latent volatility existing between the trading periods (called night volatility) has scarcely been considered in the day volatility estimation problem. The idea emerged from an observation on financial stock time series; prices at market closing differ from those at market opening the following trading day, although during the night the market is closed and thus no

transactions occur, so no intranight information exists. Despite the lack of night trades, latent (night) volatility still occurs, causing a price mismatch. We examined whether this latent night volatility can be modeled and whether, if incorporated into the conditional volatility modeling, it would help to provide better estimates of day volatility. Compared to other researchers that also modeled overnight returns, we proposed a two-factor structure in a realized GARCH setting with a GARCH equation that links day/night volatility to returns, night/day volatility, and intraday volatility of the previous day. This allowed us to retain the benefits of the Realized GARCH model of Hansen et al. (2012), namely, to take advantage of the natural relationship between the realized measure and the conditional day (and night for the models we proposed in the current paper) variance in an elegant structure that facilitates volatility estimation, allowed us to capture the return-volatility dependence, and was previously proved to outperform traditional GARCH. Below, we presented a method to capture this volatility and to insert it into the day conditional volatility equation.

The starting model is a reduced form Bivariate Realized GARCH model, which is a Realized GARCH model with night volatility information and exogenous realized measures, defined as follows:

$$r\_t = r\_t^\bullet + r\_t^\circ, \; z\_t^\circ = \frac{r\_t^\circ - \mu^\circ}{\sqrt{h\_t^\circ}}, \; z\_t^\bullet = \frac{r\_t^\bullet - \mu^\bullet}{\sqrt{h\_t^\bullet}},\tag{1}$$

$$\log h\_t^\circ = \boldsymbol{\omega}^\circ + \boldsymbol{\tau}^{(\circ 1)} (\boldsymbol{z}\_{t-1}^\bullet) + \boldsymbol{\tau}^{(\circ 2)} (\boldsymbol{z}\_{t-1}^\circ) + \boldsymbol{\beta}^\circ \cdot \log h\_{t-1}^\circ + \boldsymbol{\gamma}^\circ \cdot \log \mathbf{x}\_{t-1}.\tag{2}$$

$$\log h\_t^\bullet = \omega^\bullet + \tau^{(\bullet 1)}(z\_{t-1}^\bullet) + \tau^{(\bullet 2)}(z\_{t-1}^\circ) + \beta^\bullet \log h\_{t-1}^\bullet + \gamma^\bullet \log \mathbf{x}\_{t-1}.\tag{3}$$

where • denotes the night information, ◦ denotes the day information of the vector, *rt* is the return, *zt* <sup>∼</sup> *iid*(0, 1), *ut* <sup>∼</sup> *iid* 0, σ<sup>2</sup> *u* , *h* ◦ *<sup>t</sup>* = *var*(*r* ◦ *t* F*t*−1), *<sup>h</sup>*• *<sup>t</sup>* <sup>=</sup> *var r*• *t* |F*t*−1 F*<sup>t</sup>* = σ(*rt*, *xt*, *rt*−1, *xt*−1, ...), *r* ◦ *<sup>t</sup>* = 100 × log(*price*)*tclose* <sup>−</sup> log(*price*)*topen* , and *r*• *<sup>t</sup>* = 100 × log(*price*)*topen* <sup>−</sup> log(*price*)*t*−1*close* . As such, *rt* is the sum between night *r*• *<sup>t</sup>* and day *r* ◦ *<sup>t</sup>* returns, *z* ◦ *<sup>t</sup>* represents the standardized day returns, and *z*• *t* represents the standardized night returns, whereas μ ◦ is the means of day returns and μ• is the means of night returns. All τ's are coefficients of the standardized returns that follow to be estimated through the maximum log-likelihood function (MLE). If marked by ◦, τ represents the coefficients of the standardized returns in the equation of conditional day volatility, and if marked by •, τ represents the coefficients of the standardized returns in the equation of conditional night volatility. The numbers next to ◦ or • are for indexing purposes: For example, <sup>τ</sup>(◦1) and <sup>τ</sup>(◦2) are two coefficients of the standardized returns in the equation of conditional day volatility that follow to be estimated through MLE.

Thus, the base model is formed of three equations: The return equation, which is the sum between day (open-to-close) returns and night (close-to-open) returns, and two conditional volatility equations, as follows: The first expresses day volatility as a function of previous day (*z* ◦ *<sup>t</sup>*−1) and night (*z*• *<sup>t</sup>*−1; standardized) returns, conditional day variance (*<sup>h</sup>* ◦ *<sup>t</sup>*−1), and a realized measure of volatility (*xt*−1; realized kernel, high–low, realized variance, etc.). The second defines night volatility as a function of previous day (*z* ◦ *<sup>t</sup>*−1) and night (*z*• *<sup>t</sup>*−1; standardized) returns, conditional night variance (*h*• *<sup>t</sup>*−1), and a realized measure of volatility (*xt*−1). Notably, in this model (called reduced form for this reason), the realized measure is not endogenized nor linked to the day volatility measure through a measurement equation, but rather is treated as an exogenous variable. We added this equation to the complete form of the model that was documented in the next section. The realized measure was compounded from intraday prices recorded throughout the day.

#### *2.2. Extended Models*

We used the base model structure and extended its idea to a class of best-known GARCH-type models. We used this approach as all models used share the same structure and thus similar properties, which enabled us to set up a similar bivariate configuration. The aim was to construct a group of models that takes advantage of night volatility estimation, and also defines the existing natural relationship

between the realized measures and the conditional day and night variance. As such, we proposed four new realized models and one non-realized model: Bivariate Realized GARCH (1,1), with an endogenous component of realized measure and therefore a separate measurement equation, which we will call a complete version model; Bivariate Exponential GARCH-X (Bivariate EGARCH-X), that is a bivariate exponential generalized autoregressive conditional heteroskedastic model with an exogenous realized measure; Bivariate Realized EGARCH (1,1); Bivariate Realized GARCH (2,2); and Bivariate EGARCH (1,1). The detailed specifications of the bivariate models we propose are provided in Table 1.


**Table 1.** Summary of the bivariate realized generalized autoregressive conditional heteroskedasticity (GARCH) models proposed.

Next, we summarized the main features of each model. All share similar return equations as in the case of the base model—the daily return *rt* is the sum between open-to-close return (day return) *r* ◦ *<sup>t</sup>* and close-to-open return (night return) *r*• *<sup>t</sup>* . The GARCH equations share distinct properties but they have unique features as well. All define the day (open-to-close) volatility *h* ◦ *<sup>t</sup>* as a function of day *z* ◦ *<sup>t</sup>* and night *z*• *<sup>t</sup>* standardized returns as defined above, and also as a function of the previous day (open-to-close) volatility. Except for the Bivariate EGARCH (1,1) and the reduced form Bivariate Realized GARCH models, all other models also include the relationship between day volatility *h* ◦ *<sup>t</sup>* and intraday volatility *xt*−<sup>1</sup> in the GARCH equation. Since Bivariate EGARCH (1,1) is not a realized model, it does not contain intraday information. In our Bivariate EGARCH-X model, intraday volatility *xt*−<sup>1</sup> is treated as an exogenous variable and is thus not linked to any other variable. However, all other realized models incorporate a third equation, the measurement equation, which defines the joint dependence between

*rt* and *xt*. *xt* is thus "endogenized" by being formulated as a function of day (open-to-close) volatility, night (close-to-open) volatility, and day and night standardized returns (*z* ◦ *<sup>t</sup>* and *z*• *<sup>t</sup>* , respectively).

#### **3. Data and Estimation Methodology**

We used tick data sampled along 3537 trading days during the period of 30 August 2004–31 December 2018, corresponding to 10 stocks: AIG (American International Group, Inc.), AXP (American Express Company), BAC (Bank of America Corporation), CSCO (Cisco Systems, Inc.), F (Ford Motor Company (F)), GE (General Electric Company), INTC (Intel Corporation), JPM (JPMorgan Chase & Co.), MSFT (Microsoft Corporation), and T (AT&T Inc.). To avoid the outliers that would result from quiet days, the half trading days around the Christmas and Thanksgiving holidays were removed.

We opted for estimating intraday volatility by compounding realized kernels instead of the more widely used realized variance, as it is generally acknowledged that squared daily returns provide a poor estimation of actual intraday volatility. Realized kernels are robust for microstructure errors or frictions, which are known to cause endogenous and dependent noise terms. They are used to estimate the quadratic variation in an efficient price process when the time stamps in every day do not match (non-synchronous, with irregularly spaced observations) and when the high-frequency time series described by the prices are noisy with many microstructure effects. We compounded the realized kernels as measures of intraday volatility (*xt*) using the methodology of Barndorff-Nielsen et al. (2009, 2011). The framework is given by *Y*, a variable that is the sum of a Brownian semi-martingale and a jump process, as follows:

$$Y\_t = \int\_0^t a\_u du + \int\_0^t \sigma\_u dW\_u + f\_t. \tag{4}$$

For the purpose of our exercise, we need to find the quadratic variation of *<sup>Y</sup>*, [*Y*] = *<sup>T</sup>* <sup>0</sup> <sup>σ</sup><sup>2</sup> *udu* + *NT <sup>i</sup>*=<sup>1</sup> *<sup>C</sup>*<sup>2</sup> *<sup>i</sup>* . Barndorff-Nielsen et al. (2009, 2011) estimated it from the noisy discrete observations *X*τ*<sup>j</sup>* of *Y*τ*<sup>j</sup>* , 0 = τ<sup>0</sup> < τ<sup>1</sup> <...< τ*<sup>n</sup>* = *T*, where *X*τ*<sup>j</sup>* = *Y*τ*<sup>j</sup>* + *U*τ*<sup>j</sup>* and *U*τ*<sup>j</sup>* represents the market microstructure effects (noise). Barndorff-Nielsen et al. (2009, 2011) estimated this quadratic variation by proposing realized kernels, a non-negative estimator that is constructed as follows.

The first challenge with the tick data is the non-synchronicity. Non-synchronous trading occurs when the trades or quotes appear at irregularly spaced times across stocks, which is usually the case in stock markets, especially those with low liquidity or stale prices. Barndorff-Nielsen et al. (2011) solved this by suggesting a refresh time when all the stocks are traded. We implemented the same method by recording the prices only when (and immediately after) all of them were traded.

To eliminate start and end effects and their associated errors, which are averaged through this procedure, we proceeded to jittering (averaging) the first and last two prices, as also suggested by Barndorff-Nielsen et al. (2011)Barndorff-Nielsen et al. Having synchronized and constructed the time series by jittering at the initial and final time points, we defined the semi-definite realized kernels, as follows, according to Barndorff-Nielsen et al. (2009, 2011):

$$K(X) = \sum\_{h=-H}^{H} k(\frac{h}{H+1}) \gamma\_{h\prime} \text{ where } \gamma\_{h} = \sum\_{j=|h|+1}^{n} \mathbf{x}\_{j} \mathbf{x}\_{j-|h|} \tag{5}$$

where *k*(*x*) is a kernel weight function that has the *k*(0) = 1, *k* (0) = 0 property, and *k* is twice differentiable with continuous derivatives.

Barndorff-Nielsen et al. (2009) used a Parzen kernel as it satisfies the smoothness conditions through *k* (0) = *k* (1) = 0, and its estimates are positive. We made the same choice, and used the same Parzen kernel function:

$$k(\mathbf{x}) = \begin{cases} 1 - 6\mathbf{x}^2 + 6\mathbf{x}^3, & 0 \le \mathbf{x} \le 1/2 \\\ 2(1 - \mathbf{x})^3, & 1/2 \le \mathbf{x} \le 1 \\\ 0, \ \mathbf{x} > 1 \end{cases} \tag{6}$$

The optimal choice of bandwidth, according to Barndorff-Nielsen et al. (2009), which we chose to use, is *H*∗ = *c*∗ ξ4/5*n*3/5, with *c*<sup>∗</sup> = *k* (0) 2 *k* 0,0 • 1/5 and <sup>ξ</sup><sup>2</sup> <sup>=</sup> <sup>ω</sup><sup>2</sup> - *T T* <sup>0</sup> <sup>σ</sup><sup>4</sup> *udu* , where *c*∗ = ((12) 2 ) 1/5 <sup>=</sup> 3.5134 for

the Parzen kernel. *<sup>T</sup>* <sup>0</sup> <sup>σ</sup><sup>4</sup> *udu* is called the integrated quarticity, and, in our empirical exercise, it equals *RVsparse*. This denotes a subsampled realized variance based on 20-min returns. By calculating 1200 realized variances by shifting the first observation recorded time in 1-s increments, we obtained a number of realized variance estimators. We averaged them and obtained *RVsparse*. ω<sup>2</sup> was estimated by calculating the realized variance using every *i*th trade. We varied the starting point, and thereby produced *i* realized variances, namely *RV<sup>i</sup> dense*. Thus, our <sup>ω</sup><sup>2</sup> estimator was calculated as:

$$
\hat{\omega}\_{(j)}^2 = \frac{RV\_{\text{dense}}^{(j)}}{2n\_{(j)}}, \ j = 1, \ldots, i\_{\prime} \tag{7}
$$

where *<sup>n</sup>*(*j*) is the number of non-zero returns used to estimate *RV*(*j*) *dense*. The estimate of <sup>ω</sup><sup>2</sup> is then the average of the *j* estimates,

$$
\hat{\phi}^2 = \frac{1}{i} \sum\_{j=1}^i \hat{\phi}\_{(j)}^2. \tag{8}
$$

By design, the realized kernel is positive semi-definite and the rate of convergence is *n*1/5.

We estimated the in-sample and out-of-sample (3000th day in the sample, 24 November 2016, the cutoff point) in both the univariate and bivariate models with respect to each of the 10 stocks. The univariate models considered are the standard realized versions of the GARCH model (Realized GARCH, Realized EGARCH, EGARCH-X, and Realized GARCH (2,2)), as well as the EGARCH model. The estimated bivariate models are those mentioned in Section 2 (Bivariate EGARCH, reduced and complete forms of Bivariate Realized GARCH, Bivariate Realized EGARCH, Bivariate EGARCH-X, and Bivariate Realized GARCH (2,2)).

The estimation was performed by maximizing the total log-likelihood functions (MLE), namely the sum of partial log-likelihood functions for the returns and for the intraday measures; the ranking criterion with respect to the MLE was the partial log-likelihood function for returns solely. We used MLE to estimate both the proposed bivariate models and a number of univariate models that do not include night volatility information.

The log-likelihood function used in the estimation of the above models takes the form *l r*• *t* ,*r* ◦ *<sup>t</sup>* , *xt* <sup>=</sup> *<sup>L</sup>*<sup>1</sup> for Bivariate EGARCH and Bivariate EGARCH-X, or *l r*• *t* ,*r* ◦ *<sup>t</sup>* , *xt* <sup>=</sup> *<sup>L</sup>*<sup>1</sup> <sup>+</sup> *<sup>L</sup>*<sup>2</sup> for Bivariate Realized GARCH complete version, Bivariate Realized EGARCH (1,1), and Bivariate Realized GARCH (2,2) (Appendix A), where *L*<sup>1</sup> = −1 2 *n t*=1 ⎧ ⎪⎪⎨ ⎪⎪⎩ 2 log(2π) <sup>+</sup> log <sup>1</sup> <sup>−</sup> <sup>ρ</sup><sup>2</sup> + log *h*• *<sup>t</sup>* + log *h* ◦ *<sup>t</sup>* <sup>+</sup> (*r*• *<sup>t</sup>* − μ•) 2 /*h*• *<sup>t</sup>* +(*r* ◦ *<sup>t</sup>* − μ ◦ ) 2 /*h* ◦ *t* (<sup>1</sup> <sup>−</sup> <sup>ρ</sup>2) <sup>−</sup> <sup>2</sup><sup>ρ</sup> (1 − ρ2) (*r*• *<sup>t</sup>* <sup>−</sup> <sup>μ</sup>•)(*<sup>r</sup>* ◦ *<sup>t</sup>* − μ ◦ ) - *h*• *t h* ◦ *t* ⎫ ⎪⎪⎬ ⎪⎪⎭ and *<sup>L</sup>*<sup>2</sup> = <sup>−</sup><sup>1</sup> 2 *n t*=1 log(2π) <sup>+</sup> log σ2 *u* + *u*<sup>2</sup> *<sup>t</sup>* /σ<sup>2</sup> *u* .

To evaluate whether introducing night volatility estimations in models' equations improves the day volatility estimation, we calculated two loss functions, root mean squared error (RMSE) and mean absolute error (MAE). Based on these, we documented the number of models for each in-sample and out-of-sample estimation for each of the 10 stocks, at which MAE and RMSE were smaller. This allowed us to draw conclusions about the better performance of the bivariate or univariate models. Based on the size of the loss functions obtained at each estimation, we analyzed the performance of the new models that included night volatility estimates. This contributed to our objective by documenting whether or not night volatility information improves the estimation of day volatility with respect to the main GARCH-type of models proposed in the literature.

The maximized log-likelihood functions in univariate and bivariate estimations are provided in Tables A1 and A2 in Appendix B. As the log-likelihood functions of the bivariate models differ from those of the univariate versions (for the bivariate estimation, we maximized a bi-dimensional vector *<sup>r</sup>*• *t r* ◦ *t* with a non-null correlation factor (ρ) between its subvectors), it makes little sense to compare the values of the MLEs across the univariate and bivariate models to document an improvement or loss of performance when introducing night volatility estimates. Specifically, the log-likelihood function for the bivariate models is: log *l*(*r*• *t* ,*r* ◦ *<sup>t</sup>*) = <sup>−</sup><sup>1</sup> 2 *n t*=1 ⎧ ⎪⎪⎨ ⎪⎪⎩ 2 log(2π) <sup>+</sup> log <sup>1</sup> <sup>−</sup> <sup>ρ</sup><sup>2</sup> + log(*h*• *<sup>t</sup>* ) + log(*h* ◦ *<sup>t</sup>*) + *<sup>r</sup>*• *t* 2/*h*• *<sup>t</sup>* + *r* ◦ *t* 2/*h* ◦ *t* (<sup>1</sup> <sup>−</sup> <sup>ρ</sup>2) <sup>−</sup> <sup>2</sup><sup>ρ</sup> (1 − ρ2) *r*• *t r* ◦ *t h*• *t h* ◦ *t* ⎫ ⎪⎪⎬ ⎪⎪⎭ , where <sup>ρ</sup> <sup>=</sup> *corr r* ◦ *<sup>t</sup>* ,*r*• *t* . In the univariate models' case, the log-likelihood function is log *<sup>l</sup>*(*rt*) = <sup>−</sup><sup>1</sup> 2 *n t*=1 log(2π) <sup>+</sup> log *ht* + (*rt* <sup>−</sup> <sup>μ</sup>) 2 *ht* for EGARCH and EGARCH-X, and log *l*(*rt*) = −1 2 *n t*=1 2 log(2π) <sup>+</sup> log *ht* + (*rt* <sup>−</sup> <sup>μ</sup>) 2 *ht* <sup>+</sup> log σ2 *u* <sup>+</sup> *<sup>u</sup>*<sup>2</sup> *t* σ2 *u* for Realized EGARCH, Realized GARCH, and Realized GARCH (2,2). As such, we could not use this method to evaluate the performance of the bivariate models, as we would be comparing the values of estimations of different functions.

Thus, for the purpose of documenting the gain or loss in accuracy, we used the standard method in econometrics for evaluating the models' performance—that of calculating two loss functions (RMSE and MAE)—which would better assess whether adding night volatility information with a two-factor structure in a realized GARCH setting improves estimations of next-day volatility.

#### **4. Results**

The standard method used in econometrics to evaluate models' performance is to calculate the size of the loss functions, among which RMSE and MAE are the most common and reliable. We calculated them for both in-sample and out-of-sample estimations, and our results indicate an improvement when night volatility estimations were included in the equations of the day conditional volatility in almost every case.

We worked with a number of models that have different features and for which adding an estimation of night volatility may contribute to the volatility estimation. For example, by inspecting the results for RMSE (in-sample estimation) in Table 2, the improvement was evident for 55 out of 60 cases (1 loss function result × 6 models evaluated × 10 stocks). The cases in which the improvement could not be documented are marked with red (for RMSE) or green (for MAE) numbers in Table 2. In the five cases in which this was not evident, four of them were for Realized GARCH (2,2). This means that Realized GARCH (2,2) only shows some features that did not work better when the night volatility estimates were considered given the way in which the model was designed. This may be because, compared to the other models that model next-day volatility by only using information from the previous day and night, Realized GARCH (2,2) uses information on the previous night volatility as well as information on returns and volatility of the previous two days. We thought that this might be the problem with this model, but it would need to be proven empirically; we left this question for future work.


**Table 2.** Loss functions in univariate and bivariate estimations; in-sample.

This conclusion was strengthened by examining the MAE results. When considering MAE as an evaluation tool, the bivariate models produced superior forecasting ability in 59 out of 60 cases, indicating an improvement for the models that included night volatility estimation in the day volatility modeling. However, in only one case out of 60 was the improvement not evident, for the same Realized GARCH (2,2) model. As such, the model itself appears to be problematic, not the evaluation we performed. As mentioned above, we thought that the problem with this model was that it models conditional day volatility by including in the model information on day volatility and returns from the previous two days, instead of one day only as we did for the other models. In Bivariate Realized GARCH (2,2), we considered only one-night volatility information instead of considering the night volatility estimation from the previous two nights.

*Univ* and *Biv* stand for *Univariate* and *Bivariate*, respectively, while *com* and *red* stand for *complete* and *reduced*, respectively. Red and green numbers indicate the stances in which bivariate models perform worse than the univariate ones (when evaluated according to RMSE or MAE, respectively).

When examining the results for the out-of-sample estimations in Table 3, we found that of 60 evaluations with RMSE, 53 showed forecasting improvement when night volatility information was used. In the seven cases in which the improvement was not evident, three were recorded for the same Realized GARCH (2,2) model. The remaining four belonged to various other models, one for each. However, we observed another pattern. Most of the failures in documenting an improvement were for the same stock: AIG. This suggests that the results were sensitive not only to the model (as we explained earlier with the way in which Realized GARCH (2,2) was built), but were also sensitive to the stock choice. Since AIG persistently failed in showing an improvement when using night volatility information, AIG price recordings should be more carefully examined to understand what makes it less sensitive to this modeling suggestion, including examining the amount of the stock price differential (the difference between the market closing and the market opening prices), and also understanding the roots of the volatility transmission for this stock in particular. Again, we left this as exploratory work for the future paper. When ranked according to MAE, 58 results out of 60 indicated improvement,

whereas only two cases (among them, one for Realized GARCH (2,2)) did not. Again, both estimations indicated strong evidence in favor of including night volatility estimation in the modeling problem of day volatility.


**Table 3.** Loss functions in univariate and bivariate estimations; out-of-sample.

Counting the number of cases that fail to show improvement is valuable for two reasons: (1) It is the best tool when comparing models evaluated through MLE given that the log-likelihood functions were not similar for looking at the size of the MLE values; and (2) the cases in which we failed to see improvement indicated some consistency for a specific model and a specific stock. This opens the opportunity for future work in which we might try to understand why the Realized GARCH (2,2) model and AIG stock persistently indicated less evidence compared with other models and stocks, where by adding night volatility information, we produced improved volatility estimation.

Red and green numbers indicate the stances in which bivariate models perform worse than the univariate ones (when evaluated according to RMSE or MAE, respectively).

Thus, we concluded that the proposed bivariate models improved the forecasting performance compared with the univariate models; as such, adding night volatility estimations according to the methodology suggested improves next-day volatility estimates.

#### **5. Conclusions**

This paper provided a methodology that captures and integrates night volatility into the modeling of day volatility. In univariate context, this method led to formulating four bivariate realized GARCH models (Bivariate EGARCH-X, Bivariate Realized GARCH, Bivariate Realized GARCH (2,2), and Bivariate Realized EGARCH) and one bivariate non-realized model (Bivariate EGARCH). The novelty of this method is the incorporation of a night measure of volatility into the models, computed from price changes between the closing and opening of the trading market with a two-factor structure of the conditional variance in a realized GARCH setting that takes advantage of the natural relationship

between the realized measure and the conditional variance. This captures the leverage effect and maintains an elegant mathematical structure that facilitates the estimation of volatility.

With respect to assessing forecasting performance, the first finding was that rankings were sensitive to the stock and model choice but displayed little sensitivity to the ranking criterion and estimation methodology. However, the bivariate models were proved to perform better in most instances, compared with the univariate models. As such, we concluded that by adding night volatility estimates in the volatility models according to the methodology described, better estimates of next-day volatility could be obtained. This represents a step further from including high-frequency data in the modeling problem of the GARCH models in that estimates of night volatility are added into the equation of the day conditional variance according to the novel methodology we suggest.

The assessment to multivariate assets (e.g., portfolios of stocks) could be extended in future work by documenting a method of forecasting volatility of assets using the principal component (PC) analysis or other statistical procedures that use the orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables, taking advantage of the autoregressive conditional heteroskedastic models we proposed that use estimates of day, intraday, and night volatility. We might refer to these models as PC Bivariate Realized GARCH models and these might be used to formulate the general form of one multivariate asset's conditional variance–covariance matrix expressed in terms of conditional variances of the compounding assets and of their principal components. This would allow the estimation of the volatility of one multivariate asset through estimations of the volatility of principal components using day, intraday, and night volatility information. Then, by reducing the *n*-multivariate to a *n* − *k* stock dimension (*n* and *k* positive integers), we could estimate the new models and assess their one-day-ahead forecasting performance. Constructing models that use volatility information from the previous two days and two nights may further improve the modeling of volatility, as we noted by inspecting the results for the current bivariate form of Realized GARCH (2,2). Disseminating among the stocks according to their underlying volatility features may provide a better method of more consistently modeling their volatility patterns.

Integration of volatility estimates of highly interlinked markets that are open during the closing time of the reference market is another suggestion for further research. For example, proposing models for the U.S. market that estimate day volatility using night volatility estimates from the Asian markets open during the non-trading times of the U.S. market would allow for integration in such models of systemic risk and financial contagion related elements, with likely benefits for volatility estimation and forecasting.

**Author Contributions:** Individual contributions to the current paper were as follows: Conceptualization, M.M., X.R., and N.A.; Methodology, M.M.; Software, M.M.; Validation, M.M.; Formal Analysis, M.M.; Investigation, M.M.; Resources, M.M., X.R., and N.A.; Data Curation, M.M.; Writing—Original Draft Preparation, M.M.; Writing—Review & Editing, M.M.; Visualization, M.M.; Supervision, M.M., X.R., and N.A.; Project Administration, M.M., X.R., and N.A.; Funding Acquisition, M.M., X.R., and N.A.

**Funding:** M.M. acknowledges support from the Agency for Management of University and Research Grants (AGAUR) of the Government of Catalonia (Resolution IUE/2681/2008 of 8 August, DOGC no. 5208 of 03.09.08), ESADE Business School (Ramon Llull University), as well as from the University of Tasmania (ARC DP130100168).

**Acknowledgments:** We acknowledge the contribution of Peter Reinhard Hansen, Zhuo (Albert) Huang and Howard Howan Shek to the definition of the reduced form Bivariate Realized GARCH model described by Equations (1) to (3). We are grateful for the comments from Mardi Dungey and participants in the European Economic Association and Econometric Society Meeting 2012, Econometric Society Australasian Meeting 2013, and FIRN 2012 conference.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Disclaimer:** The views expressed herein are those of the authors and do not necessarily reflect the views of the National Bank of Romania.

#### **Appendix A. Log-Likelihood Function for the Bivariate Models**

The data are bivariate vectors compounded of two univariate vectors that refer to uncorrelated sets of information (we considered first that night volatility was uncorrelated with day volatility): *r*• *t r* ◦ *t* <sup>|</sup>*Ft*−<sup>1</sup> <sup>∼</sup> *<sup>N</sup>*(0, *<sup>h</sup>*• *<sup>t</sup>* 0 0 *h* ◦ *t* ). Accordingly, the random vector *<sup>r</sup>*• *t r* ◦ *t* depends solely on the information set available at time *<sup>t</sup>* <sup>−</sup> 1, and has a normal distribution with <sup>0</sup> 0 mean and a variance equal to the variance–covariance matrix *<sup>h</sup>*• *<sup>t</sup>* 0 0 *h* ◦ *t* . The latter is equivalent to *var r* ◦ *t* = σ ◦ *t* , *var r*• *t* = *h*• *<sup>t</sup>* and *cov r* ◦ *<sup>t</sup>* ,*r*• *t* = 0. The total volatility is given as *rt* = *r* ◦ *<sup>t</sup>* + *r*• *<sup>t</sup>* . Theory states that when a random vector (such as *<sup>r</sup>*• *t r* ◦ *t* ) is normally distributed, then its components are also normal. *<sup>r</sup>*• *t r* ◦ *t* <sup>|</sup>*Ft*−<sup>1</sup> <sup>∼</sup> *<sup>N</sup>*(0, *<sup>h</sup>*• *<sup>t</sup>* 0 0 *h* ◦ *t* ) shows that *r*• *<sup>t</sup>* |*Ft*−<sup>1</sup> , ∼ *N* 0, *h*• *t* and *r* ◦ *<sup>t</sup>* |*Ft*−<sup>1</sup> ∼ *N* 0, *h* ◦ *t* . Since a sum of two normal variables is a normal variable with the average equal to the arithmetic sum of the two component averages, *rt*|*Ft*−<sup>1</sup> ∼ *N* 0, *h*• *<sup>t</sup>* + *h* ◦ *t* , then the density function of *rt*|*Ft*−<sup>1</sup> has the form of a normal variable, that is, *f*(*rt*) = √ <sup>1</sup> σ∗ *t* √ 2π*e r*2 *t* 2*h*∗ *<sup>t</sup> ,* where *h*∗ *<sup>t</sup>* = *h*• *<sup>t</sup>* + *h* ◦ *<sup>t</sup>* is the variance of *rt*. Since n observations of *t* = 1, ... , *n* are made, the likelihood function is the ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎝ *r*1 ... *rn* ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎠ vector's density, and *r*1, ... , *rn* are independent of each other, so the likelihood function is *l*(*rt*) = \$*<sup>n</sup> <sup>t</sup>*=<sup>1</sup> *<sup>f</sup>*(*rt*) <sup>=</sup> \$*<sup>n</sup> <sup>t</sup>*=<sup>1</sup> <sup>√</sup> <sup>1</sup> *h*∗ *t* √ 2π*e r*2 *t* 2*h*∗ *<sup>t</sup>* = 1√ 2π *n* ( \$*n <sup>t</sup>*=<sup>1</sup> <sup>√</sup><sup>1</sup> *h*∗ *t* ) −1 2 *n t*=1 *r*2 *t h*∗ *<sup>t</sup>* . Taking the log of this expression and using the logarithm properties, the log-likelihood function of the total returns *rt* will become log *l*(*rt*) = log( \$*n <sup>t</sup>*=1( 1√ 2π *n* ) + log( \$*n <sup>t</sup>*=<sup>1</sup> <sup>√</sup><sup>1</sup> *h*∗ *t* ) <sup>−</sup> <sup>1</sup> 2 *n t*=1 *r*2 *t h*∗ *t* = <sup>−</sup>*<sup>n</sup>* <sup>2</sup> log(2π) <sup>−</sup> <sup>1</sup> 2 *n <sup>t</sup>*=<sup>1</sup> log *h*∗ *t* − 1 2 *n t*=1 *r*2 *t h*∗ *t* = <sup>−</sup><sup>1</sup> 2 *n t*=1 log(2π) <sup>+</sup> log *h*∗ *t* <sup>+</sup> *<sup>r</sup>*<sup>2</sup> *t h*∗ *t* . If we considered a more complete model with a non-null correlation between *r* ◦ *<sup>t</sup>* and *r*• *<sup>t</sup>* (meaning that night volatility influences day volatility), that is, *corr r* ◦ *<sup>t</sup>* ,*r*• *t* <sup>=</sup> <sup>ρ</sup> -0, the

formulation of the log-likelihood function slightly changes. Observe first that ρ does not depend on *<sup>t</sup>*, that is, the correlation is not time dependent. Then, the covariance will be *r* ◦ *<sup>t</sup>* ,*r*• *t* = *corr r* ◦ *<sup>t</sup>* ,*r*• *t* - *var r*• *t var r* ◦ *t* = ρ - *h*• *t h* ◦ *<sup>t</sup>* . This means that in the new model (with a non-null correlation), the variance–covariance matrix takes the form ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ *h*• *<sup>t</sup>* ρ - *h*• *t h* ◦ *t* ρ - *h*• *t h* ◦ *<sup>t</sup> h* ◦ *t* ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ , having the variances of *r* ◦ *t* and *r*• *<sup>t</sup>* on the first diagonal, and the covariance between *r* ◦ *<sup>t</sup>* and *r*• *<sup>t</sup>* on the second diagonal, that is *cov r* ◦ *<sup>t</sup>* ,*r*• *t* (since *cov r* ◦ *<sup>t</sup>* ,*r*• *t* <sup>=</sup> *cov r*• *t* ,*r* ◦ *t* ). As such, the <sup>|</sup>*Ft*−<sup>1</sup> conditioned distribution of the *<sup>r</sup>*• *t r* ◦ *t* vector is *<sup>r</sup>*• *t r* ◦ *t* |*Ft*−<sup>1</sup> ∼ *N*(0, ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ *h*• *<sup>t</sup>* ρ - *h*• *t h* ◦ *t* ρ - *h*• *t h* ◦ *<sup>t</sup> h* ◦ *t* ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ ). The conditional variance of *rt*, *var*(*rt*|*Ft*−<sup>1</sup> ), is *var*(*rt*|*Ft*−<sup>1</sup> ) <sup>=</sup> *var r*• *<sup>t</sup>* |*Ft*−<sup>1</sup> <sup>+</sup> *var r* ◦ *<sup>t</sup>* |*Ft*−<sup>1</sup> <sup>+</sup> <sup>2</sup>*cov r*• *<sup>t</sup>* |*Ft*−<sup>1</sup> , *r* ◦ *<sup>t</sup>* |*Ft*−<sup>1</sup> <sup>=</sup> *<sup>h</sup>*• *<sup>t</sup>* + *h* ◦ *<sup>t</sup>* + 2ρ - *h*• *t h* ◦ *<sup>t</sup>* , that is, *h*∗ *<sup>t</sup>* = *h*• *<sup>t</sup>* + *h* ◦ *<sup>t</sup>* + 2ρ - *h*• *t h* ◦ *<sup>t</sup>* . The log-likelihood function of *rt* = *r*• *<sup>t</sup>* + *r* ◦ *<sup>t</sup>* will be the same as the one iterated for the null correlation case, the only difference being that the variance encloses the correlation

term *h*∗ *<sup>t</sup>* = *h*• *<sup>t</sup>* + *h* ◦ *<sup>t</sup>* + 2ρ *h*• *t h*

◦ *t* . *Econometrics* **2019**, *7*, 41

However, we want to consider the log-likelihood function of the bivariate vector *<sup>r</sup>*• *t r* ◦ *t* and not that of the univariate vector *rt* = *r*• *<sup>t</sup>* + *r* ◦ *<sup>t</sup>* . As such, to define the new log-likelihood function, we considered the density function of the bi-dimensional normal *<sup>r</sup>*• *t r* ◦ *t* . The general form of a *p*-dimensional normal vector *Np*(μ, Σ) (a matrix with μ vector average and Σ variance–covariance matrix) takes the form *f*(*x*) = <sup>1</sup> ( √ 2π) *<sup>p</sup>* √ <sup>1</sup> *det*(Σ) *e*−1 <sup>2</sup> (*x*−μ) <sup>Σ</sup>−1(*x*−μ) , where *x* is any vector for which the density function has been calculated with p arguments, *det*(Σ) is the determinant of the variance–covariance matrix Σ, and (*x* − μ) <sup>Σ</sup>−1(*<sup>x</sup>* <sup>−</sup> <sup>μ</sup>) is the matrix product between the transpose of the (*<sup>x</sup>* <sup>−</sup> <sup>μ</sup>) vector, the inverse of matrix Σ, and the (*x* − μ) vector. As such, with *p* = 2 for the particular case of a

$$\begin{array}{ll}\text{bi-dimensional vector} \begin{pmatrix} r\_t^\bullet \\ r\_t^\circ \end{pmatrix} \text{ the density function is } f\begin{pmatrix} r\_t^\bullet \end{pmatrix} = \frac{1}{\left(\sqrt{2\pi}\right)^2} \frac{1}{\sqrt{\det(\Sigma)}} e^{-\frac{1}{2}\left(\begin{array}{c} \tau\_t^\bullet \\ r\_t^\circ \end{array}\right)} \text{ in } \Sigma\\ \begin{pmatrix} \tau\_t^\circ \end{pmatrix} \end{array}$$

which μ = 0 and Σ = ⎜⎜⎜⎜⎜⎜⎜⎝ *h*• *<sup>t</sup>* ρ *h*• *t h t* ρ - *h*• *t h* ◦ *<sup>t</sup> h* ◦ *t* ⎟⎟⎟⎟⎟⎟⎟⎠ **.** Since *det*(Σ) = *h*• *t h* ◦ *<sup>t</sup>* <sup>−</sup> <sup>ρ</sup>2*h*• *t h* ◦ *<sup>t</sup>* = *h*• *t h* ◦ *t* <sup>1</sup> <sup>−</sup> <sup>ρ</sup><sup>2</sup> , then its log ◦ 

form is log(*det*(Σ)) = log(*h*• *<sup>t</sup>* ) + log(*h <sup>t</sup>*) + log <sup>1</sup> <sup>−</sup> <sup>ρ</sup><sup>2</sup> . The inverse matrix of the variance–covariance matrix is Σ−<sup>1</sup> = <sup>1</sup> *h*• *t h* ◦ *<sup>t</sup>* (1−ρ2) ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ *h*• *<sup>t</sup>* −ρ - *h*• *t h* ◦ *t* −ρ - *h*• *t h* ◦ *<sup>t</sup> h* ◦ *t* ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ **.** As such, the product <sup>−</sup><sup>1</sup> 2 *r*• *t* ,*r* ◦ *t* Σ−<sup>1</sup> *r*• *t r* ◦ *t* becomes

−1 2 *r*• *t* ,*r* ◦ *t* Σ−<sup>1</sup> *r*• *t r* ◦ *t* = <sup>−</sup><sup>1</sup> 2 *r*• *t* 2*h* ◦ *<sup>t</sup>* +*r* ◦ *t* <sup>2</sup>*h*• *<sup>t</sup>* −2*r*• *t r* ◦ *t* ρ - *h*• *t h* ◦ *t h*• *t h* ◦ *<sup>t</sup>* (1−ρ2) . Thus, the log-likelihood function log *<sup>l</sup>*(*r*• *t* ,*r* ◦ *<sup>t</sup>*) is ◦ 

obtained by multiplying the functions *f r*• *t* ,*r t* for the *t* = 1, ... *, n*, and by taking the log of the resulting 

product log *l*(*r*• *t* ,*r* ◦ *<sup>t</sup>*) = <sup>−</sup><sup>1</sup> 2 *n t*=1 ⎧ ⎪⎪⎨ ⎪⎪⎩ 2 log(2π) <sup>+</sup> log <sup>1</sup> <sup>−</sup> <sup>ρ</sup><sup>2</sup> + log(*h*• *<sup>t</sup>* ) + log(*h* ◦ *<sup>t</sup>*) + *<sup>r</sup>*• *t* 2*h* ◦ *<sup>t</sup>* +*r* ◦ *t* <sup>2</sup>*h*• *<sup>t</sup>* −2*r*• *t r* ◦ *t* ρ *h*• *t h* ◦ *t h*• *t h* ◦ *<sup>t</sup>* (1−ρ2) ⎫ ⎪⎪⎬ ⎪⎪⎭ . By performing some simple iterations in the expression above, we obtained the final form of the bivariate log-likelihood function as log *l*(*r*• *t* ,*r* ◦ *<sup>t</sup>*) = −1 2 *n t*=1 ⎧ ⎪⎪⎨ ⎪⎪⎩ 2 log(2π) <sup>+</sup> log <sup>1</sup> <sup>−</sup> <sup>ρ</sup><sup>2</sup> + log(*h*• *<sup>t</sup>* ) + log(*h* ◦ *<sup>t</sup>*) + *<sup>r</sup>*• *t* 2/*h*• *<sup>t</sup>* +*r* ◦ *t* 2/*h* ◦ *t* (1−ρ2) <sup>−</sup> <sup>2</sup><sup>ρ</sup> (1−ρ2) *r*• *t r* ◦ *t h*• *t h* ◦ *t* ⎫ ⎪⎪⎬ ⎪⎪⎭ .

#### **Appendix B**

**Table A1.** Maximized log-likelihood functions in univariate and bivariate estimations; in-sample.



**Table A2.** Maximized log-likelihood functions in univariate and bivariate estimations; out-of-sample.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
