Hybrid CUSUM Change Point Test for Time Series with Time-Varying Volatilities Based on Support Vector Regression

This study considers the problem of detecting a change in the conditional variance of time series with time-varying volatilities based on the cumulative sum (CUSUM) of squares test using the residuals from support vector regression (SVR)-generalized autoregressive conditional heteroscedastic (GARCH) models. To compute the residuals, we first fit SVR-GARCH models with different tuning parameters utilizing a time series of training set. We then obtain the best SVR-GARCH model with the optimal tuning parameters via a time series of the validation set. Subsequently, based on the selected model, we obtain the residuals, as well as the estimates of the conditional volatility and employ these to construct the residual CUSUM of squares test. We conduct Monte Carlo simulation experiments to illustrate its validity with various linear and nonlinear GARCH models. A real data analysis with the S&P 500 index, Korea Composite Stock Price Index (KOSPI), and Korean won/U.S. dollar (KRW/USD) exchange rate datasets is provided to exhibit its scope of application.


Introduction
In this study, we aim to develop a method to detect a significant change in the conditional variance of time series with time-varying volatility using the cumulative sum (CUSUM) of squares test based on the support vector regression (SVR)-generalized autoregressive conditional heteroscedastic (GARCH) residuals. In this task, obtaining accurate predictions of volatilities is crucial because the residuals are obtained as the ratios of the observations and the forecasts of volatility. Traditionally, the prediction has been regarded as an important and challenging task in financial time series analysis owing to its non-stationary and nonlinear nature and, in particular, its volatility. The GARCH model, proposed by Engle [1] and Bollerslev [2], is the most popular model for measuring the volatility of financial time series. This model has advantages for characterizing the properties of time series well, such as persistency, heavy tails, and volatility clustering. Subsequently, several GARCH variants have been developed to better capture these properties, e.g., the asymmetric GARCH model proposed by Engle [3], the exponential GARCH (EGARCH) model proposed by Nelson [4], the threshold GARCH (TGARCH) model proposed by Zakoïan [5], the GJR-GARCH model proposed by Glosten, Jagannathan and Runkle [6], and the asymmetric power ARCH (APARCH) model proposed by Ding, Granger and Engle [7]. Refer to Carrasco and Chen [8] for the properties of these models.
Although these GARCH models perform well in general, empirical studies show that they have low forecasting performance in the presence of model misspecification. As the underlying model analysis using the S&P 500 index, KOSPI, and KRW/USD exchange rate datasets. Finally, Section 6 provides the concluding remarks.

Support Vector Regression for the GARCH Model
Let us consider the GARCH model: where g is an unknown function, p, q are nonnegative integers, and { t } is a sequence of independent and identically distributed (i.i.d.) random variables with zero mean and unit variance. When g is linear, the inferences for stationary GARCH(p, q) models are well developed in the literature (Francq and Zakoïan [39] and Francq and Zakoïan [40]). Ideally, this approach provides an excellent analytic tool for prediction, but the accuracy of prediction can be deteriorated owing to a violation of the assumptions such as stationarity and linearity. As such, in this study, we use the SVR-GARCH method for prediction. SVR is a revision of the SVM, initially proposed by Cortes and Vapnik [41], which seeks an optimal hyperplane that separates the inputs by maximizing the margins between the support vectors and the hyperplane. Equipped with a high generalization ability (Abe [42]), a notable strength of the SVM is that its usage can be extended to nonlinear prediction and classification methods because of its ability to map the inputs to a higher dimensional feature space via nonlinear kernel functions (Cortes and Vapnik [41]).
SVR aims to identify a nonlinear function of the form: where x denotes a vector of inputs, w and b are vectors of the regression parameter, and φ is a known kernel function. In the context of this paper, x is comprised of squared inputs up to lag p (y 2 t−1 , . . . , y 2 t−p ) and conditional variance up to lag q (σ 2 t−1 , . . . , σ 2 t−q ). Notice that w and b have the same dimension as x. The optimal w and b that yield the best approximation of f are determined by exploiting the structure of the -insensitive loss function (Vapnik [17]): which neglects the errors lying inside the -band surrounding the estimated function f . The SVR aims to seek a function that keeps the data inside such a band and dismisses excessive complexity (Smola and Schölkopf [18]). Given the input vectors x i , scalar output y i , i = 1, . . . , n, and a constant C > 0, we construct the objective function of the SVR using the -insensitive loss function (Vapnik [17]): subject to where ξ 1,i , ξ 2,i > 0 denote slack variables that allow some points to lie outside the -band with a penalty and C denotes a trade-off between the function complexity and the training error. Notice that the former two constraints designate the upper and lower bound, respectively.
To obtain the optimal w and b, we formulate an unconstrained optimization problem using Lagrange multipliers (Smola and Schölkopf [18]): Then, the optimal solution must satisfy the following Karush-Kuhn-Tucker conditions (Vapnik [17]): and: Substituting Equations (5) and (6) into (4), we obtain the following dual form: subject to where α 1,i and α 2,i denote dual variables (Vapnik [17]). Consequently, the optimization problem in (7) yields the solutionsŵ,b,f of w, b, f of the following form (Smola and Schölkopf [18]): . Note thatŵ is determined uniquely, whereasb is not, although the two cases rarely coincide. One way to determineb is to obtain the average of the above two values (Abe [42]). Furthermore, when constructing a model with the SVR here, we employ the Gaussian kernel for K in (8): As such, we have to determine the tuning parameters γ 2 , C in (3), and in the loss function (2) appropriately. Accordingly, we consider a cube of (C, γ 2 , ) with 1 ≤ C ≤ 100, 0.1 ≤ γ 2 ≤ 1, and 0.1 ≤ ≤ 1. We then employ a grid search method on this cube.
An analogous approach can be adopted for the GARCH time series (Chen, Härdle and Jeong [13]) by constructing the design matrix X with each row x s comprising the lagged time series of lag p and the lagged conditional variance of lag q. Specifically, for the GARCH(p, q) model, To reflect the structure above, we present a five-step procedure to obtain the predicted SVR-GARCH function.
Given the time series of length n + n and the space of tuning parameters, the procedure is illustrated as follows:

•
Step 1. Prescribe the points to be evaluated within this space, then divide the given time series into training and validation time series of size n and n , respectively. This preliminary procedure is required for the subsequent task of validating the fitted SVR-GARCH model, which determines the best tuning parameter sets.

•
Step 2. Note that the conditional variance σ 2 t of (9) is unknown. As a remedy, replace σ 2 t with the initial estimatesσ 2 t , which plays the role of a proxy of σ 2 t . The estimateσ 2 t is based on the training time series using a moving average method (Niemira [43]): where m is a positive integer. When t is smaller than m,σ 2 t is computed as an average of the first to the tth squares of observations, i.e.,σ 2 Step 3. Given a set of tuning parameters, we estimate g in (1) withĝ using the SVR with σ 2 t replaced byσ 2 t . Then, the estimateσ 2 t of σ 2 t is obtained as: Step 4. Applying the estimated SVR-GARCH model and using the same proxy formula as in Step 2 for the validation time series, the mean absolute error (MAE) is computed as follows: The MAE escalates the robustness of the model against outliers and therefore provides more flexibility in a model fitting than the root mean squared error.

•
Step 5. Repeat Steps 2 to 4 for all the tuning parameter sets selected in Step 1 and choose the combination that minimizes the MAE. Then, perform Steps 2 and 3 using the training and validation time series together to determine the final model, which is used in obtaining the residuals.

Remark 1.
One might contemplate utilizing different methods regarding the selection of tuning parameters, such as the random search method proposed by Bergstra and Bengio [44]. However, we employ the grid search method in our simulation study because it outperforms the random search method. Furthermore, the choice of m in the construction of proxy volatilities is critical in constructing a stable test. In this study, we report the case of m = 5 because it provides reasonably satisfactory results the most consistently.

Remark 2.
One may consider iteratively updatingσ 2 t until a specific convergent criterion is satisfied, as observed in Chen, Härdle and Jeong [13] and Lee, Lee and Moon [38]. However, theĝ obtained from this iterative approach makes a function either "flatten out" after each iteration, ultimately yielding a constant function, or drastically follow the peak outliers among the initialσ 2 t 's. Thus, this approach is inadequate and disregarded in our study.

Remark 3.
Instead of MAE, other loss functions, such as the mean squared error (MSE), the root mean squared error (RMSE), and the mean absolute percentage error (MAPE), could be employed. These loss functions, however, are likely to yield suboptimal results. To illustrate, when either the MSE or the RMSE is used, the result deteriorates as these amplify or shrink the losses according to their values. On the other hand, MAPE yields inconsistent results because it uses unobservableσ 2 t rather thanσ 2 t .

Hybrid CUSUM Test via the SVR-GARCH Model
Based on the work of Lee, Tokutsu and Maekawa [35], Oh and Lee [34] developed the CUSUM of squares test for the GARCH(1,1) model as follows: where ω > 0, α ≥ 0, β ≥ 0, α + β < 1, and t are i.i.d. random variables with zero mean and unit variance, to perform a test for a parameter change in ω, α, β. More precisely, the null and alternative hypotheses are formulated as follows: H 0 : ω, α, β remain the same over t = 1, . . . , n.vs.
To test these, ω, α, β are estimated using their consistent estimatorsω,α,β, such as Gaussian quasi-maximum likelihood estimates (QMLEs) in Francq and Zakoïan (2004), and the residualsˆ t are recursively obtained via the equation: with the initials y 0 = 0, σ 2 0 = y 1 . Subsequently, the CUSUM of squares test is given bŷ Letting with τ 2 = Var( 2 1 ), Oh and Lee [34] verified that under regularity conditions,T n behaves asymptotically the same as T n under H 0 as n tends to ∞, and thus, the limiting null distribution ofT n is the same as sup 0≤s≤1 |B • (s)|, where B • denotes a Brownian bridge on the unit interval because T n converges to the supremum of a Brownian bridge in distribution due to Donsker's invariance principle; see Billingsley [45]. The critical values then can be obtained asymptotically. For instance, the null hypothesis H 0 of no changes is rejected ifT n ≥ 1.3397 at the level of 0.05, which can be obtained through Monte Carlo simulations. Furthermore, provided that a change point exists, the location of change is identified as:k Multiple change points can be detected by following the scheme of Inclán and Tiao [24]. The same approach can be adopted for more general nonlinear location-scale time series models as seen in Oh and Lee [32].
This CUSUM procedure can be applied to the residuals obtained through the SVR-GARCH models, as can be seen in our simulation study, since they will quite likely behave like the true errors t . Lee, Lee and Moon [38] recently adopted the SVR-ARMA scheme to calculate the residuals in the change point detection problem on ARMA models and affirmed its validity.
In the next section, we evaluate the SVR-GARCH model in the previous section for various settings, wherein we only consider the case of p = q = 1 as most financial time series can be sufficiently described with GARCH(1,1) models (Hansen and Lunde [46]).

Simulation Results
In this section, we evaluate the performance of the SVR-GARCH model through simulations using the linear, asymmetric, threshold, GJR, and exponential GARCH models. For this task, we generate a time series of length n = 500, 1000, and 2000 to evaluate the empirical sizes and powers at the level of 0.05. The experiment is comprised of the following four steps.

•
Step 1. Generate a time series of length 2n from a prescribed GARCH model.

•
Step 2. Follow the estimation scheme described in Section 3 with m = 5. In this procedure, the first 0.7n number of time series constitute the training set, and the following 0.3n number of time series constitute the validation set.

•
Step 3. Conduct the CUSUM of squares test described in Section 3. We utilize the remaining n number of time series as a testing set.

•
Step 4. Repeat Steps 1 to 3 1000 times iteratively, and then, compute the empirical sizes and powers.
In the power evaluation, we assume that the change point existed in the middle of the time series. Here, we utilize R 3.5.1 running on Windows 10 and the packages "e1071" and "fGarch".
When examining empirical power, we first consider the change in parameters. To illustrate this, we consider the change in ω, the sum of α and β for the GARCH, TGARCH, and log-linear GARCH models, and the change in b for the AGARCH model. Similarly, for the GJR-GARCH model, we consider a change in the sum of α 1 , α 2 , and β instead. Under the null hypothesis, the parameters are set as follows: 1.
GARCH model: AGARCH model: TGARCH model: log-linear GARCH model: Tables 1 and 2 indicate the results for the GARCH and AGARCH models, respectively. In both models, no size distortions are observed in most cases. Furthermore, the results confirm that the test performs well in general, in terms of power. Furthermore, as anticipated, the power increases as the sample size increases. Tables 3-5 depict the results for the GJR-GARCH, TGARCH, and log-linear GARCH models, respectively. Although some mild distortions of size can be noticed in the GJR-GARCH model, the results are mostly similar to those of the GARCH and AGARCH models.
For the log-linear GARCH model, the truncated residuals are employed in the construction of the CUSUM test, that is,˜ t :=ˆ t I(|ˆ t | < N) + N I(|ˆ t | ≥ N) with N = 20 to mitigate the influence of extreme outliers. This modification does not hamper the asymptotic behavior of the CUSUM of squares test, whereas it significantly improves the performance of the CUSUM test in terms of stability and power, as demonstrated in Table 5.  Table 5. Empirical sizes and powers for the log-linear GARCH(1,1) model. Finally, to simulate the situation in which the underlying model is unknown, we examine the case where both the parameters and the underlying model itself change. Analogously to the simulation study regarding the above log-linear GARCH model, we truncate the residuals and implement˜ t to formulate the CUSUM statistic in Cases 1 and 2: 1.
AGARCH(1,1) changes to GJR-GARCH(1,1). Table 6 shows that the test produces powers comparable to those with only parameter changes. Overall, our findings strongly support the validity of the SVR-GARCH model for detecting a change. Table 6. Empirical powers for model changes (log-linear GARCH = log-GARCH).

Real Data Analysis
In this section, we analyze the log-returns of the daily index of the S&P 500 from 3 January 2012 to 30 September 2016, the Korea Composite Stock Price Index (KOSPI) from 2 October 2012 to 28 June 2019, and the exchange rate between U.S. dollars (USD) and South Korean Won (KRW), denoted as KRW/USD, from 2 January 2012 to 30 September 2016. We obtain these time series from the website "investing.com". The two South Korean economic indices are selected because they are well appreciated to be susceptible to various international affairs of the country owing to its geographical location and heavy dependence on the exports of the South Korean economy. Prior to fitting the SVR-GARCH model, we first inspect the autocorrelation function (ACF) and partial autocorrelation function (PACF) of the log-returns of each training time series, consisting of the first half of the entire time series, to examine the presence of irregular patterns of autocorrelations. This step is needed to verify the adequacy of the training time series. Figure 1 shows that the ACFs and PACFs for all three training datasets support stationarity to a great extent, indicating that estimation via the SVR-GARCH model with these training time series would not undermine the outcomes. The selected training datasets are highlighted with yellow shaded areas in   Figure 2, the original index, shows a relatively consistent trend of steady increase before the change point. However, a double-dip is observed after the change point, which coincides with the fact that a change in trend is often accompanied by a change in variation in stock markets. Furthermore, Figures 3 and 4 plot the result of the change point detection process for the KOSPI and KRW/USD indices, respectively. Here, we obtain the CUSUM statistic values of 1.5081 and 4.2754, respectively, which indicates the detection of a change in both cases at the level of 0.05. The identified location of a change point for the KOSPI index appears to be 22 January 2018, whereas that of KRW/USD is 26 September 2014. Similar to the case of the S&P 500 index, a significant change in trend is also observed in the original datasets, as shown in the plot (a) of Figures 3 and 4. Our finding also illuminates that the log-return of the KRW/USD index experiences a more significant change in volatility, in contrast to the other two cases.
In general, it is not feasible to find the matching incidents that cause change points in economic indices either because they can be obscured from media or they may affect the market with a significant delay. Nevertheless, for the KOSPI example, we could deduce that unstable international affairs could be the culprit behind such changes. These affairs include the nuclear weapon test of North Korea in September 2017, and the participation of North Korea in the Pyeongchang Winter Olympic Games in February 2018, which caused massive turmoil in the Korean peninsula. However, one could argue that the raising of the interest rate by the FRB twice in March and June 2018, was a much more significant factor affecting both the changes in volatility and trend, considering that the impact of the international affairs of the Korean peninsula has been limited on many occasions. In contrast, in the case of the KRW/USD index, the report from the Bank of Korea in January 2015, reasoned that the high volatility after the change point was due to a byproduct of ending quantitative easing and the accompanying improvement of the U.S. economy.

Conclusions
In this study, we proposed the CUSUM of squares test based on the residuals obtained with the SVR-GARCH model in order to detect a parameter change in the volatility of time series. Monte Carlo simulations were conducted with various linear and nonlinear GARCH models, including the GARCH, GJR-GARCH, AGARCH, TGARCH, and log-linear GARCH models, and the obtained results confirmed the validity of the SVR-GARCH method. Our method was then applied to the analysis of financial datasets such as the S&P 500, KOSPI, and KRW/USD indices and detected one change in all cases. Overall, our findings supported the validity of our method and the practicality in financial time series analysis. Here, we only considered a plain SVR method for emphasizing the hybrid of the CUSUM and SVR methods and for easy access to general readers. However, more sophisticated methods could be employed for refinement concerning the selection of tuning parameters and kernel functions, although they do not necessarily guarantee a better performance, possibly due to over-fitting; see Remark 2 of Lee, Lee and Moon [38]. Furthermore, in this study, we only investigated the retrospective change point test. However, a similar method could be applied to an on-line monitoring process (Huh, Oh and Lee [47]), which aims at an early detection of an anomaly in sequentially observed time series. In this case, the training data could be redesigned via a rolling window procedure. As these issues go beyond the scope of the current study, we leave them to our future project.