Forecasting Realized Volatility Using a Nonnegative Semiparametric Model

This paper introduces a parsimonious and yet flexible semiparametric model to forecast financial volatility. The new model extends the linear nonnegative autoregressive model of Barndorff-Nielsen and Shephard (2001) and Nielsen and Shephard (2003) by way of a power transformation. It is semiparametric in the sense that the distributional and functional form of its error component is partially unspecified. The statistical properties of the model are discussed and a novel estimation method is proposed. Simulation studies validate the new method and suggest that it works reasonably well in finite samples. The out-of-sample forecasting performance of the proposed model is evaluated against a number of standard models, using data on S&amp;P 500 monthly realized volatilities. Some commonly used loss functions are employed to evaluate the predictive accuracy of the alternative models. It is found that the new model generally generates highly competitive forecasts.


Introduction
Financial market volatility is an important input for asset allocation, investment, derivative pricing, and financial market regulation. Not surprisingly, how to model and forecast financial volatility has been a subject of extensive research. Numerous survey papers are now available on the subject, with hundreds of reviewed research articles. Excellent survey articles on the subject include Bollerslev, Chou, and Kroner (1992), Bollerslev, Engle, and Nelson (1994), Ghysels, Harvey, and Renault (1996), Poon and Granger (2003), and Shephard (2005).
In this vast literature, ARCH and stochastic volatility (SV) models are popular parametric tools. These two classes of models are motivated by the fact that volatilities are time-varying.
Moreover, they offer ways to estimate past volatility and forecast future volatility from return data. In recent years, however, many researchers have argued that one could measure latent volatility by realized volatility (RV), see e.g. Andersen, Bollerslev, Diebold, and Labys (2001) (ABDL hereafter) and Barndorff-Nielsen and Shephard (2002), and then build a time series model for volatility forecasting using observed RV, see e.g. Andersen, Bollerslev, Diebold, and Labys (2003). An advantage of this approach is that "models built for the realized volatility produce forecasts superior to those obtained from less direct methods" (ABDL, 2003). In an important study, ABDL (2003) introduced a new Gaussian time series model for logarithmic RV (log-RV) and established its superiority for RV forecasting over some standard methods based on squared returns. Their choice of modeling log-RV rather than raw RV is motivated by the fact that the logarithm of RV, in contrast to RV itself, is approximately normally distributed. Moreover, conditional heteroskedasticity is greatly reduced in log-RV.
Following this line of thought, in this paper we introduce a new time series model for RV.
For the S&P 500 monthly RV, we show that although the distribution of log-RV is closer to a normal distribution than that of raw RV, normality is still rejected at all standard significance levels. Moreover, although conditional heteroskedasticity is reduced in log-RV, there is still evidence of remaining conditional heteroskedasticity. These two limitations associated with the logarithmic transformation motivate us to consider a more flexible transformation, that is, the so-called Tukey's power transformation which is closely related to the well-known Box-Cox transformation. In contrast to the logarithmic transformation, Tukey's power transformation or the Box-Cox transformation is generally not compatible with a normal error distribution as the support for the normal distribution covers the entire real line. 1 This well-known truncation problem further motivates us to use nonnegative error distributions. The new model, which we call a Tukey nonnegative type autoregression (TNTAR), is flexible, parsimonious, and has a simple forecast expression. Moreover, the numerical estimation of the model is very fast and can easily be implemented using standard computational software.
The new model is closely related to the linear nonnegative models described in Barndorff- Nielsen and Shephard (2001) and Nielsen and Shephard (2003). In particular, it generalizes the discrete time version of the nonnegative Ornstein-Uhlenbeck process of Barndorff-Nielsen and Shephard (2001) by (1) applying a power transformation to volatility; (2) leaving the dependency structure and the distribution of the nonnegative error term unspecified. Our work is also related to Yu, Yang, and Zhang (2006) and Gonçalves and Meddahi (2011) where the Box-Cox transformation is applied to stochastic volatility and RV, respectively. The main difference between our model specification and theirs is that an unspecified (marginal) distribution with nonnegative support, instead of the normal distribution, is induced by the transformation. Moreover, our model is loosely related to Higgins and Bera (1992), Hentschel (1995) and Duan (1997) where the Box-Cox transformation is applied to ARCH volatility, and to Fernandes and Grammig (2006) and Chen and Deo (2004). Finally, our model is related to a recent study by Cipollini, Engle, and Gallo (2006) where an alternative model with nonnegative errors is used for RV. The main difference here is that the dynamic structure for the transformed RV is linear in our model, whereas the dynamic structure for the RV is nonlinear in theirs.
Our proposed model is estimated using a two-stage estimation method. In the first stage, a nonlinear least squares procedure is applied to a nonstandard objective function. In the second stage a linear programming estimator is applied. The finite sample performance of the proposed estimation method is studied via simulations.
The TNTAR model is used to model and forecast the S&P 500 monthly RV and its out-ofsample performance is compared to a number of standard time series models previously used in the literature, including the exponential smoothing method and two logarithmic long-memory ARFIMA models. Under various loss functions, we find that our parsimonious nonnegative model generally generates highly competitive forecasts. While this paper considers the application of forecasting RV, there are a number of applications beyond financial data for which our model may be useful. For example, modeling and forecasting climatological or telecommunication time series may be interesting alternative applications for our nonnegative model. While our model is related to several models in the literature, to the best of our knowledge, our specification is new in two ways. First, it is based on Tukey's power transformation. Second, the distribution and functional form of its error component are partially unspecified. Moreover, the estimation method that we propose is new.
The rest of the paper is organized as follows. Section 2 motivates and presents the new model.
In Section 3 a novel estimation method is proposed to estimate the parameters of the new model.
In Section 4 the finite sample performance of the new method is studied via simulations. Section 5 describes the S&P 500 realized volatility data and the empirical results. In the same section we also outline the alternative models for volatility forecasting and present the loss functions used to assess their forecast performances. Finally, Section 6 concludes.

A Nonnegative Semiparametric Model
Before introducing the new TNTAR model, we first review two related time series models previously used in the volatility literature, namely, a simple nonnegative autoregressive (AR) model and the Box-Cox AR model. Barndorff-Nielsen and Shephard (2001) introduced the following continuous time model for fi-

Related volatility models
In the above z is a Lévy process with independent nonnegative increments, which ensures the positivity of σ 2 (t) (see Equation (2) in Barndorff-Nielsen and Shephard, 2001). Applying the Euler approximation to the continuous time model in (1) yields the following discrete time model where ϕ = 1 − λ and u t+1 = z(λ(t + 1)) − z(λt) is a sequence of independent identically distributed (i.i.d.) random variables whose distribution has a nonnegative support. A well known nonnegative random variable is the generalized inverse Gaussian, whose tails can be quite fat.
Barndorff- Nielsen and Shephard (2001) discuss the analytical tractability of this model. In the case when u t+1 is exponentially distributed, Nielsen and Shephard (2003) derive the finite sample distribution of a linear programming estimator for ϕ for the stationary, unit root and explosive cases. 2 Simulated paths from model (2) typically match actual realized volatility data quite well.
See e.g. Figure 1(c) in Barndorff- Nielsen and Shephard (2001). Unfortunately, so far little empirical evidence establishing the usefulness of this model has been reported.
Two restrictions seem to apply to model (2). First, since its errors are independent, conditional heteroskedasticity is not allowed for. The second restriction concerns the ratio of two successive volatilities. More specifically, from (2) it can be seen that σ 2 t+1 /σ 2 t is bounded from below by ϕ, almost surely, implying that σ 2 t+1 cannot decrease by more than 100(1 − ϕ)% compared to σ 2 t . Since the AR parameter ϕ of the model typically is estimated using linear programming, in practice, this restriction is automatically satisfied. For instance, the full sample estimate of ϕ in our empirical study is 0.262, implying that σ 2 t cannot decrease by more than 73.8% from one time period to the next. Indeed, 73.8% is the maximum percentage drop in successive monthly volatilities in the sample, which took place on November 1987.
In a discrete time framework, a popular parametric time series model for volatility is the 2 See Section 3 for a detailed discussion on the linear programming estimator. lognormal SV model of Taylor (2007) given by where r t is the return, σ 2 t is the latent volatility, and ε t and t are two independent Gaussian noises. In this specification volatility clustering is modeled as an AR(1) for the log-volatility. The logarithmic transformation in (4) serves three important purposes: First, it ensures the positivity of σ 2 t . Second, it removes heteroskedasticity. Third, it induces normality. Yu et al. (2006) introduced a closely related SV model by replacing the logarithmic transformation in Taylor's volatility equation (4) with the more general Box-Cox transformation (Box and Cox, 1964), where Compared to the logarithmic transformation, the Box-Cox transformation provides a more flexible way to improve normality and reduce heteroskedasticity. A nice feature of the Box-Cox AR model given by (5)-(6) is that it includes several standard specifications as special cases, including the logarithmic transformation (λ = 0) and a linear specification (λ = 1). In the context of SV, Yu et al. (2006) find empirical evidence against the logarithmic transformation. Chen and Deo (2004) and Gonçalves and Meddahi (2011) are interested in the optimal power transformation. In the context of RV, Gonçalves and Meddahi (2011) find evidence of non-optimality for the logarithmic transformation. They further report evidence of negative values of λ as the optimal choice for various data generating processes. Our empirical results reinforce this important conclusion, although our approach is very different.
While the above discrete time models have proven useful for modeling volatility, there is little documentation on their usefulness for forecasting volatility. Moreover, the Box-Cox transformation is known to be incompatible with a normal error distribution. This is the well-known truncation problem associated with the Box-Cox transformation in the context of Gaussianity.

Realized volatility
In the ARCH or SV models, volatilities are estimated parametrically from returns observed at the same frequency. In recent years, however, it has been argued that one can measure volatility in a model-free framework using an empirical measure of the quadratic variation of the underlying efficient price process, that is, RV. RV has several advantages over ARCH and SV models. First, by treating volatility as directly observable, RV overcomes the well known curse-of-dimensionality problem in the multivariate ARCH or SV models. Second, compared to the squared return, RV provides a more reliable estimate of integrated volatility. This improvement in estimation naturally leads to gains in volatility forecasting.
Let RV t denote the RV at a lower frequency (say daily or monthly) and p(t, k) denote the log-price at a higher frequency (say intra-day or daily). Then RV t is defined by where N is the number of higher frequency observations in a lower frequency period. 3 The theoretical justification for RV as a volatility measure comes from standard stochastic process theory, according to which the empirical quadratic variation converges to integrated variance as the infill sampling frequency tends to zero (ABDL, 2001, Barndorff-Nielsen and Shephard, 2002, Jacod, 2017. The empirical method inspired by this consistency has recently become more popular with the availability of high-frequency data. In a recent important contribution, ABDL (2003) find that a Gaussian long-memory model for the logarithmic daily realized variance provides more accurate forecasts than the GARCH(1,1) model and the RiskMetrics method of J.P. Morgan (1996). The logarithmic transformation is used since it is found that the distribution of logarithmic realized variance, but not raw realized variance, is approximately normal. In Table 1 we report (to 3 decimals) some summary statistics for monthly RV, log-RV and power-RV for the S&P 500 data in our empirical study over the period Jan 1946-Dec 2004, including the skewness, kurtosis, and p-value of the Jarque-Bera test statistic for normality. 4 For RV, the departure from normality is overwhelming. While the distribution of log-RV is much closer to a normal distribution than that of RV, there is still strong evidence against normality.
To compare the conditional heteroskedasticities, in Figure 1 we plot squared OLS residuals ( 2 it , i = 1, 2, 3), obtained from AR(1) regressions for RV, log-RV and power-RV, respectively, against each corresponding explanatory variable (lagged RV, log-RV and power-RV). For ease of comparison, superimposed are smooth curves fitted using the LOESS method. It is clear that while the logarithmic transformation reduces the conditional heteroskedasticity there is still evidence of 3 In ABDL (2003) RV is referred to as the realized variance, ∑ N k=2 [p(t, k) − p(t, k − 1)] 2 . Although the authors build time series models for the realized variance, they forecast the realized volatility. In contrast, the present paper builds time series models for, and forecasts, the realized volatility, which seems more appropriate. Consequently, the bias correction, as described in ABDL (2003), is not required. 4 The power parameter is -0.278 which is the estimate of λ in our proposed TNTAR model obtained using the entire S&P 500 monthly RV sample. See sections 3 and 5 for further details. it in the residuals. The power transformation further reduces the conditional heteroskedasticity of RV. While the logarithmic transformation reduces the impact of large observations (extreme deviations from the mean), the second plot of Figure 1 suggests that it is not as effective as anticipated. In contrast, the power transformation with a negative power parameter is able to reduce the impact of large observations further. Thus, the results indicate that there is room for further improvements over the logarithmic transformation. A more detailed analysis of the S&P 500 data is provided in Section 5.

The model
In this paper, our focus is on modeling and forecasting RV. To this end, let us first consider the RV version of model (5), where t is a sequence of independent N(0, σ 2 ) distributed random variables and h(x, λ) is given by (6).
If λ = 0, we may rewrite (8) as where RV λ t is a simple power transformation. A special case of (9) is a linear Gaussian AR(1) model, obtained when λ = 1: If λ = 0 in (8), we have the log-linear Gaussian AR(1) model previously used in the literature: While the specification in (8) is more general than the log-linear Gaussian AR(1) model (11), it has a serious drawback. In general, solving for RV t , the right hand side of (9) has to be nonnegative with probability one or almost surely (a.s.). This requirement is violated since a normal error distribution has a support covering the entire real line.
This drawback motivates us to explore an alternative model specification for RV. Our proposed nonnegative TNTAR model is of the form with the power parameter λ = 0, AR parameter ϕ > 0, and (a.s.) positive initial value RV 1 . The errors u t driving the model are nonnegative, possibly non-i.i.d., random variables. In the simplest case, u t is assumed to be a sequence of m-dependent, identically distributed, continuous random variables with nonnegative support [η, ∞), for some unknown η ≥ 0. 5 It is assumed that m ∈ N is finite and potentially unknown. Hence, the distribution and functional form of u t is partially unspecified. We expect ϕRV λ t−1 to be the dominating component in (12) and do not model u t parametrically.
The power transformation RV λ t is closely related to John W. Tukey's ladder of power transformations for linearizing data (Tukey, 1977), partially illustrated in (13) below: The nonnegative restriction on the support of the error distribution ensures the positivity of RV λ t . Hence, our model does not suffer from the truncation problem of the classical Box-Cox model (8). As the distribution of u t is left unspecified, some very flexible tail behavior is allowed for. Consequently, the drawback in the Box-Cox AR model (8) is addressed in the proposed TNTAR model (12).
In the classical Box-Cox model, the transformation parameter λ is required to induce linearity and normality and at the same time eliminate conditional heteroskedasticity. These are too many requirements for a single parameter. In our model, the role of the Tukey-type power transformation is to improve linearity and reduce conditional heteroskedasticity, not to induce normality. To illustrate this, suppose that a square root transformation is applied with λ = 1/2 in (12), then RV t = ϕ 2 RV t−1 + 2ϕ √ RV t−1 u t + u 2 t and the conditional variance of raw RV is time-varying. 6 An intercept in the model is superfluous because the support parameter η can be strictly positive.
Our model echoes (8), with the normal distribution replaced by a nonnegative distribution. If λ = 1 and its errors are i.i.d., our model becomes the discrete time version of Equation (2) in Barndorff-Nielsen and Shephard (2001). In general, the distributional and functional form is not assumed to be known for the error component. Hence, the TNTAR model combines a parametric component for the persistence with a nonparametric component for the error. On the one hand, the new model is highly parsimonious. In particular, there are only two parameters that need to be estimated for the purpose of volatility forecasting, namely ϕ and λ. On the other hand, the specification is sufficiently flexible for modeling the error.
As mentioned earlier, there exists a lower bound for the percentage change in volatility in model (2). A similar bound applies to our model. It is easy to show that RV t /RV t−1 ≤ ϕ 1/λ if λ < 0 (upper bound), and RV t /RV t−1 ≥ ϕ 1/λ if λ > 0 (lower bound). Typical estimated values of ϕ and λ in (12) for our empirical study are 0.639 and -0.278, respectively, suggesting that RV t cannot increase by more than 500% from one time period to the next. As we will see later, our proposed estimator for λ depends on the ratios of successive RV's and hence the bound is endogenously determined.
6 More generally, suppose that λ = 1/n for some natural number n, then RV t = ϕ n RV t−1 + u t n = ∑ n k=0 ( n k )ϕ n−k RV In this section we consider the estimation of the parameters ϕ and λ, and a one-step-ahead forecast expression, for the TNTAR model. First, we consider the special case when λ is assumed to be known. Some common power transformations include λ = 1/n (the nth root transformation) and its reciprocal, λ = −1/n. Second, we consider the more general case when both ϕ and λ are unknown and need to be estimated. We then examine the finite sample performance of the proposed estimation method via simulations.

Robust estimation of ϕ
If the true value of the power transformation parameter is known, a natural estimator for ϕ in (12) given the sample RV 1 , . . . , RV T of size T and the nonnegativity of the errors is The estimator ϕ T in (14) can be viewed as the solution to a linear programming problem. Because of this, we will refer to it as a linear programming estimator (LPE). This estimator is also the conditional (on RV 1 ) maximum likelihood estimator (MLE) of ϕ when the errors in (12) exponentially distributed random variables, cf. Nielsen and Shephard (2003). Interestingly, the LPE is strongly consistent for more general error specifications, including heteroskedasticity and m-dependence. It is robust in the sense that its consistency conditions allow for certain model misspecifications in u t . For example, the order of m-dependence in the error sequence and the conditional distribution of RV t may be incorrectly specified. Moreover, the LPE is strongly consistent even under quite general forms of heteroskedasticity and structural breaks. For a more detailed account of the properties of the LPE, see Preve (2015).
Like the ordinary least squares (OLS) estimator for ϕ, the LPE is distribution-free in the sense that its consistency does not rely on a particular distributional assumption for the error component. However, the LPE is in many ways superior to the OLS estimator. For example, its rate of convergence can be faster than O p (T −1/2 ) even for ϕ < 1, whereas the rate of covergence for the OLS estimator is faster than O p (T −1/2 ) only for ϕ ≥ 1, see Phillips (1987). Furthermore, unlike the OLS estimator the consistency conditions for the LPE do not involve the existence of any higher order moments.
Under additional technical conditions, Davis and McCormick (1989) and Feigin and Resnick (1992) obtain the limiting distribution of a LPE for which (14) appear as a special case when λ = 1 and the errors are i.i.d.. The authors show that the accuracy of the LPE depends on the index of regular variation at zero (or infinity) of the error distribution function. For example, for standard exponential errors, the index of regular variation at zero is 1 and the LPE converges to ϕ at the rate of O p (T −1 ). In general, a difficulty in the application of the limiting distribution is that the index of regular variation at zero appears both in a normalizing constant and in the limit. Datta and McCormick (1995) avoid this difficulty by establishing the asymptotic validity of a bootstrap scheme based on the LPE.
It is readily verified that the LPE in (14) is positively biased and stochastically decreasing in T, that is, ϕ < ϕ T 2 ≤ ϕ T 1 a.s. for any T 1 < T 2 . 7 Hence, the accuracy of the LPE either remains the same or improves as the sample size increases (cf. Figure 2).
To illustrate the robustness of the LPE, consider a covariance stationary AR(1) model under the possible misspecification where t is a sequence of non-zero mean i.i.d. random variables. For m > 0 the (identically distributed) errors u t are serially correlated. In this setting the OLS estimator for ϕ is inconsistent while the LPE remains consistent. In the first panel of Figure 2 we plot 100 observations simulated from the nonnegative ARMA(1,1) model, RV t = ϕRV t−1 + t + ψ t−1 with ϕ = 0.5, ψ = 0.75 and standard exponential noise. In the second panel of Figure 2 we plot the sample paths of the recursive LPEs and OLS estimates for ϕ obtained from the simulated data. In each iteration, a new observation is added to the sample used for estimation. It can be seen that the LPEs quickly approach the true value ϕ, whereas the OLS estimates do not. Moreover, the OLS estimates fluctuate much more than the LPEs when the sample size is small, suggesting that the LPE is less sensitive to extreme deviations from the mean than the OLS estimator in small samples.
We now list simple assumptions under which the consistency of the LPE in (14) holds. More general assumptions, allowing for an unknown number of unknown breaks in the error mean and variance, under which the LPE converges to ϕ for a known λ are given in Preve (2015).
Assumption 1. The power transformation parameter λ = 0 in (12) is known. The AR parameter ϕ > 0, and the initial value RV 1 is a.s. positive. The errors u t driving the autoregression form a sequence of m-dependent, identically distributed, nonnegative continuous random variables. The order, m, of the dependence is finite.
Assumption 1 allows for various kinds of m-dependent error specifications, with m ∈ N potentially unknown. For example, serially correlated finite-order MA specifications. Since the functional form and distribution of u t are taken to be unknown, the formulation is nonparametric.
It is important to point out that Assumption 2 is satisfied for any error distribution with unbounded nonnegative support. 7 Whenever necessary we use the subscript T to emphasize on the sample size. Theorem 1. Suppose that assumptions 1-2 hold. Then the LPE in (14) is strongly consistent for ϕ in (12). That is, ϕ T converges to ϕ a.s. as T tends to infinity.
The convergence of ϕ T is almost surely (and, hence, also in probability). Our interest is in forecasting raw RV, not the power transformation of RV in (12). Let RV T+1 denote a forecast of RV T+1 made at time T. A simple approximation to the optimal mean squared error, one-stepahead, forecast of RV T+1 at time T is given by the sample average converges to u i in distribution as T tends to infinity under assumptions 1-2.

Estimation of ϕ and λ
In practice, we usually do not know the true value of λ. In this section we propose an LPE based two-stage estimation method for ϕ and λ in the TNTAR model (12). In doing so, we also establish a general expression for its one-step-ahead forecast. The estimators are easily computable using standard computational software such as Matlab.
Joint estimation of ϕ and λ is non-trivial, even under certain parametric and simplifying assumptions for u t . For example, even in the simple case when u t is a sequence of independent exponentially distributed random variables it appears that the MLEs of ϕ and λ are inconsistent.
Because of this we propose an estimation method based on the LPE for ϕ.
In our LPE based two-stage estimation method, we first choose λ T to minimize the sum of squared one-step-ahead prediction errors: where respectively. Although our estimator for λ looks like the standard nonlinear least squares (NLS) estimator of Jennrich (1969), the two approaches are quite different because in our model an explicit expression for E(RV t | RV t−1 ) is not available. In fact, the NLS estimators of λ and ϕ, that minimizes ∑ T t=2 (RV l t − pRV l t−1 ) 2 , always take values of 0 and 1, respectively, and hence are inconsistent.
The intuition behind the proposed estimation method is that we expect RV t ( λ T ) to be close to E(RV t | RV t−1 ) for large values of T. This is not surprising since the TNTAR model (12) implies In the second stage, we use the LPE to estimate ϕ. More specifically, While we minimize the sum of squared one-step-ahead prediction errors when estimating λ, other criteria, such as minimizing the sum of absolute one-step-ahead prediction errors, can be used. We have experimented with absolute prediction errors using the S&P 500 data and found that our out-of-sample forecasting results for the TNTAR model are quite insensitive to the choice of the objective function in the estimation stage. However, the objective function with squared prediction errors performs better in simulations.
It is beyond the scope of this paper to derive asymptotic properties for the two-stage estimators. However, under primitive assumptions, the consistency of λ T and ϕ T can be established using the fundamental consistency result for extremum estimators. Moreover, under high-level assumptions, the martingale central limit theorem can be used to establish the asymptotic distribution of λ T .
With an estimated λ and ϕ, a general one-step-ahead semiparametric forecast expression for the TNTAR model is given by Of course, in line with Granger and Newbold (1976), several forecasts of RV T+1 may be considered. For example, one could base a forecast on the well known approximation E[h(y)] ≈ h[E(y)] using h(y) = y 1/λ . However, this approximation does not take into account the nonlinearity of h(y). 8

Monte Carlo Studies
We now examine the performance of our estimation method via simulations. We consider two experiments in which data are generated by the nonnegative model with i.i.d. standard exponential driving noise t .
In the first Monte Carlo experiment λ is assumed to be known and we only estimate ϕ using the LPE in (14). In this case the consistency is robust to the first-order moving average specification of u t . Hence, we simulate data from the model with the value of ψ being different from zero.  Table 3: Simulation results for the proposed two-stage estimation method. Summary statistics for λ T and ϕ T based on data generated by the nonnegative process RV λ t = ϕRV λ t−1 + t with i.i.d. standard exponential noise t . Bias and MSE denotes the empirical bias and mean squared error, respectively. Results based on 100, 000 Monte Carlo replications. second experiment. Here the negative bias arises because λ is estimated. In sum, it seems that the proposed estimation method works well, especially when T is reasonably large.

An Empirical Study
We also study the performance of the proposed model for forecasting actual RV relative to popular existing models. Before we report empirical results, we first review some alternative models and criteria to evaluate the performance of different models.

Alternative models
Numerous models and methods have been applied to forecast stock market volatility. For example, ARCH-type models are popular in academic publications and RiskMetrics is widely used in practice. Both methods use returns to forecast volatility at the same frequency. However, since the squared return is a noisy estimator of volatility ABDL (2003) instead consider RV and present strong evidence to support time series models based directly on RV in terms of forecast accuracy. We also compare the performance of our model against the exponential smoothing method, a RV version of RiskMetrics. The AR and log-AR models are defined by (10) and (11), respectively. We now review the exponential smoothing method, the ARFIMA model, and the HAR model.

Exponential smoothing
Exponential smoothing (ES) is a simple method of forecasting, where the one-step-ahead forecast of RV T+1 at time T is given by with 0 < α < 1. The exponential smoothing formula can be understood as the RV version of RiskMetrics, where the squared return, r 2 T , is replaced by RV T . Under the assumption of conditional normality of the return distribution, r 2 T is an unbiased estimator of σ 2 t . RiskMetrics recommends α = 0.94 for daily data and α = 0.97 for monthly data.
To see why the squared return is a noisy estimator of volatility even under the assumption of conditional normality of the return distribution, suppose that r t follows (3). Conditional on σ t , it is easy to show that (Lopez, 2001) This implies that with a probability close to 0.74 the squared return is at least 50% greater, or at most 50% smaller, than the true volatility. Not surprisingly, Andersen and Bollerslev (1998) find that RiskMetrics is dominated by models based directly on RV. For this reason, we do not use RiskMetrics directly. Instead, we use (17) with α = 0.97, which assigns a weight of 3% to the most recently observed RV. We remark that the forecasting results of Section 5 were qualitatively left unchanged when other values for α were used.

ARFIMA(p, d, q)
Long range dependence is a well documented stylized fact for volatility of many financial time series. Fractional integration has previously been used to model the long range dependence in volatility and log-volatility. The autoregressive fractionally integrated moving average (ARFIMA) was considered as a model for logarithmic RV in ABDL (2003) and Deo, Hurvich, and Lu (2006), among others. In this paper, we consider two parsimonious ARFIMA models for log-RV, namely, an ARFIMA(0, d, 0) and an ARFIMA(1, d, 0).
The ARFIMA(p, d, 0) model for log-RV is defined by where the parameters µ, β 1 , . . . , β p and the memory parameter d are real valued, and ε t is a sequence of independent N(0, σ 2 ε ) distributed random variables. Following a suggestion of a referee, we estimate all the parameters of the ARFIMA model using an approximate ML method by minimizing the sum of squared one-step-ahead prediction errors. See Beran (1995), Chung and Baillie (1993), and Doornik and Ooms (2004) for detailed discussions about the method and for Monte Carlo evidence supporting it. Compared to the exact ML method of Sowell (1992), there are two advantages to the approximate ML method. First, it does not require d to be less than 0.5. Second, it has smaller finite sample bias. Compared to the semi-parametric methods, it is also more efficient. 9 The one-step-ahead forecast of RV T+1 at time T of an ARFIMA(p, d, 0) for log-RV with p = 0 is given by and with p = 1 by and Γ(·) denotes the gamma function.

HAR
The HAR model proposed by Corsi (2009) is one of the most popular models for forecasting volatility. Given that we will forecast monthly RV in the empirical study, we modify the original HAR model with monthly, quarterly and yearly components. The original HAR model was proposed to model daily RV. We apply the modified model to raw RV (HAR) and to log-RV (log-HAR). The model for raw RV can be expressed as where the parameters β 0 , . . . , β 3 are real valued, RV t is the realized volatility of month t, and RV m t−1 = RV t−1 , RV q t−1 = 1 3 ∑ 3 i=1 RV t−i , RV y t−1 = 1 12 ∑ 12 i=1 RV t−i denote the monthly, quarterly and 9 We also applied the exact ML method of Sowell (1992) and the exact local Whittle estimator of Shimotsu and Phillips (2005) in our empirical study and found that the forecasts remained essentially unchanged. yearly lagged RV components, respectively. This specification of RV parsimoniously captures the high persistence observed in our empirical study. The one-step-ahead forecast of RV T+1 at time T is given by The corresponding forecast of the HAR model in (19) for log-RV is where σ 2 is the estimated variance of the independent N(0, σ 2 ) distributed errors t .

Forecast accuracy measures
It is not obvious which accuracy measure is more appropriate for the evaluation of the outof-sample performance of alternative time series models. Rather than making a single choice, we use four measures to evaluate forecast accuracy, namely, the mean absolute error (MAE), the mean absolute percentage error (MAPE), the mean square error (MSE), and the mean square percentage error (MSPE). Let RV it denote the one-step-ahead forecast of RV t at time t − 1 of model i and define the accompanying forecast error by e it = RV t − RV it . The four accuracy measures are defined, respectively, by where P is the length of the forecast evaluation period.
An advantage of using MAE instead of MSE is that it has the same scale as the data. The MAPE and the MSPE are scale independent measures. For a comprehensive survey on these and other forecast accuracy measures see Hyndman and Koehler (2006).
When calculating the forecast error, it is implicitly assumed that RV t is the true volatility at time t. However, in reality the volatility proxy RV t is different from the true, latent, volatility.
Several recent papers discuss the implications of using noisy volatility proxies when comparing volatility forecasts under certain loss functions. See, for example, Andersen and Bollerslev (1998), Hansen and Lunde (2006) and Patton (2011). The impact is found to be particularly large when the squared return is used as a proxy for the true volatility, but diminishes with the approximation error. In this paper, the true (monthly) volatility is approximated by the RV using 22 (daily) squared returns. As a result, the approximation error is expected to be considerably smaller than in the case of using a single squared return.

Data
The data used in this paper consists of daily closing prices for the S&P 500 index over the period January 2, 1946-December 31, 2004, covering 708 months and 15,054 trading days. We measure the monthly volatility using realized volatility calculated from daily data. Denote the log-closing price on the k'th trading day in month t by p(t, k). Assuming there are T t trading days in month t, we define the monthly RV as where 1/T t serves the purpose of standardization.
In order to compare the out-of-sample predictive accuracy of the alternative models, we split the time series of monthly RV into two subsamples. The first time period is used for the initial estimation. The second period is the hold-back sample used for forecast evaluation. When computing the forecasts we use a recursive scheme, where the size of the sample used for parameter estimation successively increases as new forecasts are made. The time series plot of monthly RV for the entire sample is shown in Figure 3, where the vertical dashed line indicates the end of the initial sample period used for estimation in our first forecasting exercise.  Table 4: Summary statistics for the S&P 500 monthly RV data. JB is the p-value of the Jarque-Bera test under the null hypothesis that the data are from a normal distribution,ρ i is the ith sample autocorrelation.
kurtosis is 28.791 indicating that the distribution of RV is non-Gaussian. In contrast, log-RV has a much smaller kurtosis (3.657) and is less skewed (0.389). It is for this reason that we include Gaussian time series models for log-RV in the exercise. However, a formal test for normality via the JB statistic rejects the null hypothesis of normality of log-RV, suggesting that further improvements over log-linear Gaussian approaches are possible.
Higher order sample autocorrelations are in general slowly decreasing and not statistically negligible, indicating that RV and log-RV are predictable. To test for possible unit roots, augmented Dickey-Fuller (ADF) test statistics were calculated. The ADF statistic for the sample from 1946 to 2004 is -5.69 for RV and -5.43 for log-RV, which is smaller than -2.57, the critical value at the 10% significance level. Hence, we reject the null hypothesis that RV or log-RV has a unit root.

Empirical results
Each alternative model was fitted to the in-sample RV data and used to generate one-step-ahead out-of-sample forecasts. 10 Following a suggestion of a referee, we also included a standard GARCH(1,1) (sGARCH) and a realized GARCH(1,1) with a log-linear specification (realGARCH), Hansen, Huang, and Shek (2012). 11 Since a forecast frequency of one month is sufficiently important in practical applications, we focus on one-step-ahead forecasts in this paper. However, multi-step-ahead forecasts can be obtained in a similar manner.
We perform two out-of-sample forecasting exercises. In both exercises, we use the recursive scheme, where the size of the sample used to estimate the alternative models grows as we make forecasts for successive observations. 12 More precisely, in the first exercise, we first estimate all the alternative models with data from the period January 1946-June 1975 and use the estimated models to forecast the RV of July 1975. We then estimate all models with data from January 10 The Ox language of Doornik (2009) was used to estimate the two ARFIMA models. Matlab code and data used in this paper can be downloaded from http://www.mysmu.edu/faculty/yujun/research.html.
11 The sGARCH and realGARCH models were estimated using monthly log-returns and the rugarch R package of Ghalanos (2019).
12 While we consider the recursive forecasting scheme one could, of course, also consider the rolling or fixed scheme.

Sample including the 1987 crash
In the first exercise, the first month for which an out-of-sample volatility forecast is obtained is July 1975. In total 354 monthly volatilities are forecasted, including the volatility of October 1987 when the stock market crashed and the RV is 0.026.
In Figure 4, we plot the monthly RV and the corresponding one-month-ahead TNTAR forecasts for the out-of-sample period, July 1975 to December 2004. It seems that the TNTAR model captures the overall movements in RV reasonably well. The numerical computation of the 354 forecasts is fast and takes less than five minutes on a standard desktop computer.
In Figure 5, we plot the recursive estimates, λ T and ϕ T . While λ T takes values from -0.45 to -0.28, ϕ T ranges between 0.58 and 0.64. It may be surprising to see that the path of ϕ T is non-monotonic. This is because the estimates of the power transformation parameter, λ, are varying over time. Our empirical estimates of λ seem to corroborate well with the optimal value of λ obtained by Gonçalves and Meddahi (2011) using simulations in the context of a GARCH diffusion and a two factor SV model. While ϕ T is quite stable, λ T jumps in October 1987.
For comparison, we also consider a TNTAR model with λ taken to be known. Visual inspection, see Figure 6, shows that a power transformation with λ = −1/2 improves linearity considerably. 13 We denote the corresponding TNTAR model TNTAR * , and employ the LPE based forecasting scheme proposed in Preve (2015): We first fit the TNTAR model using the LPE and calculate LP residuals Due to the robustness of the LPE, simple semiparametric forecasts in the (possible) presence of structural breaks are then obtained by applying a one-sided moving median. More specifically, as a simple one-month-ahead forecast we take RV T+1 = m T , where m T is the sample median of , the reciprocals of the by ϕ T / √ RV T shifted, squared last 12 LP residuals.

Sample post the 1987 crash
To examine the sensitivity of our results with respect to the 1987 crash and the 1997 crash due to the Asian financial crisis, we redo the forecasting exercise so that the first month for which an out-of-sample volatility forecast is obtained is January 1988 and the last month is September 1997.  In Figure 7, we plot the monthly RV and the corresponding one-month-ahead TNTAR forecasts for the out-of-sample period, January 1988-September 1997. As before, forecasts from the TNTAR model captures the overall movements in RV reasonably well.

Concluding Remarks
In this paper, a simple time series model is introduced to model and forecast RV. The new TNTAR model combines a nonnegative valued process for the error term with the flexibility of Tukey's power transformation. The transformation is used to improve linearity and reduce heteroskedasticity while the nonnegative support of the error distribution overcomes the truncation problem in the classical Box-Cox setup. The model is semiparametric as the order of m-dependence, support parameter η and functional form of its error term are left unspecified. Consequently, the proposed model is highly parsimonious, having only two parameters that need to be estimated for the purpose of forecasting. A two-stage estimation method is proposed to estimate the parameters of the new model. Simulation studies validate the new estimation method and suggest that it works reasonably well in finite samples.
We empirically examine the forecasting performance of the proposed model relative to a number of existing models, using monthly S&P 500 RV data. The out-of-sample performances were evaluated under four different forecast accuracy measures (MAE, MAPE, MSE and MSPE).
We found empirical evidence that our nonnegative model generates highly competitive volatility forecasts.
Why does the simple nonnegative model generate such competitive forecasts? Firstly, as shown in Section 2.2, the logarithmic transformation may not reduce heteroskedasticity and improve normality as well as anticipated. A more general transformation may be required. Secondly, the nonnegative model is highly parsimonious. This new approach is in sharp contrast to the traditional approach which aims to find a model that removes all the dynamics in the original data. When the dynamics are complex, a model with a rich parametrization is called for. This approach may come with the cost of over-fitting and hence may not necessarily lead to superior forecasts. By combining a parametric component for the persistence and a nonparametric error component, our approach presents an effective utilization of more recent information.
Although we only examine the performance of the proposed model for predicting S&P 500 realized volatility one month ahead, the technique itself is quite general and can be applied in many other contexts. First, the method requires no modification when applied to intra-day data to forecast daily RV. In this context, it would be interesting to compare our method to the preferred method in ABDL (2003). Second, our model can easily be extended into a multivariate context by constructing a nonnegative vector autoregressive model. Third, while we focus on stock market volatility in this paper, other financial assets and financial volatility from other financial markets can be treated in the same fashion. Fourth, as two alternative nonnegative models, it would be interesting to compare the performance of our model with that of Cipollini et al. (2006). Finally, it would be interesting to examine the usefulness of the proposed model for multi-step-ahead forecasting. These extensions will be considered in later work.