==== Front Entropy (Basel) Entropy (Basel) entropy Entropy 1099-4300 MDPI 33287077 10.3390/e22111312 entropy-22-01312 Article Monitoring Volatility Change for Time Series Based on Support Vector Regression https://orcid.org/0000-0003-1109-6768Lee Sangyeol * Kim Chang Kyeom Kim Dongwuk Department of Statistics, Seoul National University, Seoul 08826, Korea; ckdk95@snu.ac.kr (C.K.K.); bitsaem@snu.ac.kr (D.K.) * Correspondence: sylee@stats.snu.ac.kr; Tel.: +82-2-880-8814 17 11 2020 11 2020 22 11 131221 10 2020 13 11 2020 © 2020 by the authors.2020Licensee 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/).This paper considers monitoring an anomaly from sequentially observed time series with heteroscedastic conditional volatilities based on the cumulative sum (CUSUM) method combined with support vector regression (SVR). The proposed online monitoring process is designed to detect a significant change in volatility of financial time series. The tuning parameters are optimally chosen using particle swarm optimization (PSO). We conduct Monte Carlo simulation experiments to illustrate the validity of the proposed method. A real data analysis with the S&P 500 index, Korea Composite Stock Price Index (KOSPI), and the stock price of Microsoft Corporation is presented to demonstrate the versatility of our model. GARCH-type time seriesCUSUM monitoringsupport vector regressionparticle swarm optimization ==== Body 1. Introduction In this paper, we study the cumulative sum (CUSUM) monitoring procedure to sequentially detect a significant change in time series with conditional volatilities. Since [1,2], the CUSUM method has been an acclaimed tool to detect an anomaly among observations in statistical process control (SPC). In SPC, a control chart is a primary component that graphically describes the behavior of sequentially observed time series. For examples of SPC, see [3] for control charts, [4] for the analysis of complex environmental time series, and [5] for backcasting-forecasting time series. In particular, the CUSUM chart is one of the most frequently adopted methods among various fields of SPC. For an overview of SPC, we refer to [6]. In general, the performance of control charts is measured with average run length (ARL). However, instead of the conventional control charts, some authors, such as [7], alternatively took the approach of controlling type I errors in probability instead of controlling ARL to deal with the monitoring process in autoregressive time series. This design of the sequential monitoring method has merit in its ability to attain a lower false alarm rate, as seen in [8], who took a similar approach to dealing with generalized autoregressive conditional heteroscedastic (GARCH) time series. For more references as to the monitoring process in time series, see [9]. Here, inspired by the previous studies, we aim to hybridize the CUSUM monitoring scheme with support vector regression (SVR) for GARCH-type time series. Support vector machine (SVM) is one of the most popular nonparametric learning methods used predominantly for classification and regression [10]. In particular, compared with traditional model-based methods, the support vector regression (SVR) is more advantageous to approximating the nonlinearity of the underlying dynamic structure of datasets. Refer to [11,12,13,14,15], and the papers cited therein. Moreover, SVR is structured to exploit the quadratic programming optimization problem, and to implement the structural risk minimization principle [16]. This facilitates SVR to be an estimation scheme that optimally minimizes the empirical risk while being relatively parsimonious [17]. This further enables the SVR to perform adequately for variously sized datasets, including the case where its size is relatively small. Refs. [18,19] recently used the SVR with a hybridization of the CUSUM method to detect a change point in SVR-autoregressive and moving average (ARMA) and SVR-GARCH models. Therein, the residual-based CUSUM test has been adopted by referring to the previous studies of [20,21,22]. See [23] for a general overview of the change point detection problem using CUSUM methods. In the same spirit, here we also take the approach of the CUSUM of squares test, rather than the score vector-based CUSUM test used in [7,8], as the former largely outperforms the latter in terms of stability and power, as seen in the empirical study of [21]; above all, the latter is only available for the model-based CUSUM test. As pointed out in the precedent works of [18,19], the role of the accurately computed residuals is substantially important, and thereby, so is the precise prediction of the used SVR method. As the optimal choice of tuning parameters matters to a significant degree as well, we consider adopting particle swarm optimization (PSO) [24] for the purpose of tuning parameter optimization. For a theoretical background, we refer to [25,26,27]. See also the survey papers of [28,29,30]. These are population-based algorithms and have been widely used to obtain optimal tuning parameters in SVR. Therefore, we evaluate the performance of the proposed SVR-GARCH-based CUSUM monitoring method equipped with the PSO through Monte Carlo experiments. The rest of this paper is organized as follows. Section 2 introduces the CUSUM monitoring and explains its fundamental principle and application to GARCH-type time series. Section 3 describes SVR and PSO in general, then elaborates on the monitoring process for SVR-GARCH models in more detail. Section 4 presents Monte Carlo simulations conducted to evaluate the performance of the proposed method. Section 5 performs a real data analysis using the S&P 500 index, KOSPI, and the stock returns of Microsoft Corporation datasets. Finally, Section 6 provides concluding remarks. 2. CUSUM Monitoring Procedure In this section, we introduce our monitoring process starting from the independent and identically distributed (iid) sample case. Let us consider the problem of monitoring an anomaly in the variance from a stream of observations ϵt of mean zero up to time n. Under the null hypothesis of no anomalies, we assume that ϵt are iid with unit variance over time t=1,…,n. Namely, we test H0:Var(ϵt)=1,t=1,…,nvs.H1:Var(ϵt)=1,t=1,…,kVar(ϵt)≠1,t=k+1,…,nforsomek=2,…,n−1. For this task, we first define (1) Wk=∑t=1k(ϵt2−1)/τ, for each k≥1, where τ2=Var(ϵ12), which is assumed to be known. Ref. [7] considered the monitoring process based on: (2) Tn(1)=max1≤k≤nTn(k)=max1≤k≤nmaxm≤k1nWm−Wk, and later [8] additionally considered the monitoring process based on: (3) Tn′=max1≤k≤nTn′(k)=max1≤k≤nsup1≤m0, at some point k0=2,…,n−1, |minm≤kWm−Wk| for k>k0 would tend to have larger values due to the shift, which, however, is not true for maxm≤kWm−Wk, indicating that Tn(2) would be able to detect the shift well while Tn(1) would not do so; by contrast, if  Eϵk2 has a negative shift, Tn(1) would be able to detect the shift well while Tn(2) would not do so. Additionally, Tnmax is preferable over Tn′ in (3) as the latter tends to detect well a change that occurs in the middle of time series. Using Donsker’s invariance principle [31] and the fact that sup0≤s≤tB(s)−B(t)=|B(t)| in distribution for any standard Brownian motion B, as described in [7], we have that as n→∞, Tn(1)→DT=sup0≤t≤1(sup0≤s≤tB(s)−B(t))=Dsup0≤t≤1|B(t)|,Tn(2)→DT′=sup0≤t≤1(sup0≤s≤t(−B)(s)−(−B)(t))=Dsup0≤t≤1|B(t)|, where →D denotes a convergence in distribution, and X=DY implies the distribution functions of random variables X and Y being equivalent. Furthermore, by the continuous mapping theorem, Tnmax in (4) satisfies (6) Tnmax→DT∨T′. Then, the null hypothesis is rejected if Tnmax is larger than a constant c, which is determined asymptotically as the number satisfying P(T∨T′>c)=α for any significance α∈(0,1). In practice, we can obtain the critical value c empirically using Monte Carlo simulations. For example, c=2.46509 when α=0.05. Thus, an anomaly is signaled at k when Tnmax(k)=Tn(1)(k)∨Tn(2)(k)>c for some k=1,…,n. This monitoring process can be extended to the GARCH(1,1) model [32]: (7) yt=σtϵt,t∈Z,σt2=ω+αyt−12+βσt−12 with ω>0,α≥0,β≥0 satisfying α+β<1 and Eϵ14<∞. This model is widely used to model financial time series with high volatilities as it captures well the volatility clustering phenomenon. For its properties and characteristics and various GARCH variants, see [19,33], and the papers cited therein. Here, the monitoring process can be constructed based on residuals ϵ^t=yt/σ^t0 with σ^t02=ω0+α0yt−1,02+β0σ^t−12 with some initial values y0 and σ0 when the true parameters ω0,α0,β0 are known in advance from past experience. However, when they are unknown, which is a prevalent phenomenon in practice, we instead construct the CUSUM test with residuals ϵ^t=yt/σ^t with σ^t2=ω^+α^yt−12+β^σ^t−12, where ω^,α^,β^ are estimators of ω,α,β obtained from a given training sample. Then, we employ the following: (8) T^nmax=T^n(1)∨T^n(2), where T^n(i) are the same as Tn(i) in (2) and (5), i=1,2, with Wk in (1) replaced by W^k=∑t=1k(ϵ^t2−ϵ^¯2)/τ^, and ϵ^¯2 and τ^2 are the sample mean and variance of the residuals obtained from the training sample, respectively. We reject the null hypothesis if T^nmax>c for the aforementioned critical value c>0. Theoretically, the limiting distribution of T^nmax is anticipated to be the same as the one in (6) when certain conditions are fulfilled, namely m=m(n) satisfies m/n→∞ as n→∞ (e.g., m∼cnloglogn for some c>0) and (ω^,α^,β^) is a (m-consistent) Gaussian quasi-maximum likelihood estimator (QMLE) as in [33], the proof of which is rather standard and omitted for brevity. The monitoring process of the nonlinear time series via SVR is provided in Section 3.3 below. 3. Monitoring Procedure via SVR-GARCH Model 3.1. Support Vector Regression Support vector regression (SVR) is a nonparametric function estimating method, which is a branch of the support vector machine that originated from [34] as a classification tool. Here, SVR is utilized to conduct CUSUM tests, as it effectively incorporates circumstances where the underlying structure of the conditional variance presented in (7) is unknown, and possibly nonlinear. Its details will be elaborated in Section 3.3 below. In SVR, we seek to find a function of the form: f(x)=〈w,ϕ(x)〉+b, where x denotes a vector of explanatory variables, w and b are regression parameters to be estimated, and ϕ is an implicit kernel operator that satisfies K(x1,x2)=ϕ(x1)ϕ(x2) for some pre-defined kernel function K. To obtain the estimates of the regression parameters, we then formulate the problem where we can exploit its structure to employ a quadratic programming method [17]: (9) minimize12||w||2+C∑i=1n(ξi+ξi*), subjecttoyi−{wTϕ(xi)+b}≤ϵ+ξi*{wTϕ(xi)+b}−yi≤ϵ+ξiξi≥0,ξi*≥0, where ξi,ξi*>0 denote slack variables, C>0 denotes a penalty term, and ϵ is a tuning parameter that determines the level of tolerance regarding the error. Note that C and ϵ are user-defined tuning parameters of the model. By the duality of the Karush–Khun–Tucker condition, we obtain the following optimization problem with respect to the Lagrange multipliers αi and αi* as dual variables [16]: (10) maximize−12∑i=1n∑j=1n(αi−αi*)(αj−αj*)ϕ(xi)Tϕ(xj)−ϵ∑i=1n(αi+αi*)+∑i=1n(αi−αi*)yi, subjectto∑i=1n(αi−αi*)=00≤αi≤C,0≤αi*≤C. Consequently, the solutions of the optimization problem (10) constructs the estimate f^ of f as follows: w^=∑i=1n(αi−αi*)ϕ(xi),f^(x)=∑i=1n(αi−αi*)K(xi,x)+b^, where b^ can be obtained with various approaches. See [17] or [35] for references. The tuning parameters of our SVR-GARCH model consist of C,ϵ, and γ2, as we utilize the Gaussian kernel function K(x,y)=exp−||x−y||22γ2. Moreover, as the tuning parameters are user-defined, it must either be selected prior to employing the SVR, or be optimized. Here, we utilize the PSO to select the set of tuning parameters. The procedure of the PSO algorithm is elucidated in more detail in the next section. 3.2. Particle Swarm Optimization PSO is one of the widely acclaimed meta-heuristic methodologies and is extensively adopted when optimizing tuning parameters in various machine learning techniques. Initially proposed by [24], its ability of optimization in nonlinear and nonconvex problems is inspired by the movement of organisms in bird flocks or animal herds. Here, we refer to [30] to describe the procedure. Given a suitable objective function, a standard PSO algorithm aims to optimize a d-dimensional parameter x in the following search space: X={x=(x1,…,xd)T∈Rd:lk≤xk≤uk,k=1,⋯d},d≥1, for some uk,lk∈R, k=1,⋯d. A particle j at time t, denoted by a d-dimensional vector xj(t) in X with xj(t)=(xj1(t),xj2(t),⋯,xjd(t))T∈X, is considered as a candidate of the set of optimal solutions. Moreover, a swarm at time t is defined as S(t)=(x1(t),⋯,xN(t))T. Each particle j has its own d-dimensional velocity vector at time t, denoted by vj(t)∈V for j=1,⋯N, where V={v=(v1,…,vd)T∈Rd:vmin≤vk≤vmax,k=1,⋯,d} is a velocity space with vmin and vmax being the lower and upper bound of each velocity element, respectively [36]. Furthermore, the best optimal solution pj(t) is defined as: pj(t)=(pi1,⋯,pjd)T∈X(j,t), where X(j,t):={xj(1),⋯,xj(t)} denotes the trajectory of the particle j until time t. Then, the global best optimal solution of the trajectory of all particles until time t is defined as g(t)=(g1,⋯,gd)T∈{pj(1),…,pj(t)}. In each generation at time t∈1,Tmax, where Tmax is a prescribed maximum generation time, the j-th particle in swarm S(t) is updated as follows: vj(t+1)=w(t)vk(t)+c1z1(pi(t)−xi(t))+c2z2(g(t)−xi(t))xi(t+1)=xi(t)+vi(t+1), where c1,c2 are acceleration factors in R and z1,z2 are random variables generated from a uniform distribution U[0,1]. In addition, w(t) is an inertia term at time t, which is calculated as follows: w(t)=wstart−wendTmax−tTmax+wend, where wstart and wend are predefined lower and upper bounds of the inertia values, respectively. Ref. [30] recently proved that the optimal solution obtained from the PSO algorithm converges to the global optimum with probability 1 under regularity conditions. The process of the aforementioned PSO algorithm is condensed in Algorithm 1 below. Algorithm 1 Standard PSO algorithm 1: procedure PSO(N,wmax,wmin,c1,c2,Tmax) 2: Acceleratedfactor=(c1,c2),maxgeneration=Tmax 3:  Initializelocationandvelocity; 4:  while t0 is some initial value. This alteration improves both the stability and the performance of the model, as portrayed in Section 4. Moreover, to further enhance the stability of the estimation process of g^, we log-transform the response variable of (11), then take the exponential to obtain σ^t2. To elaborate, we recursively obtain σ^t2 through γt:=logσt2 as follows: (12) γt=g*(yt−12,σt−12),g*=logg,γ^t=g^*(yt−12,σ˜t−12)←γ˜t=logσ˜t2,σ^t2=exp(γ^t). Given g^ and a space of tuning parameters Θ, we then employ the PSO algorithm to obtain an optimal set of tuning parameters θ*∈Θ by evaluating the mean absolute error (MAE): MAE=1m−l∑t=l+1m|σ^t2−σ˜t2|, where σ^t2 and σ˜t2 are obtained with a validation time series. Ultimately, the finalized g^ is obtained by utilizing all the training samples {y1,…,ym} and θ*. The test set of length n emerges sequentially in practice, denoted by ym+1,…,ym+n. Then, utilizing g^, we yield the residuals ϵ^t=yt/σ^t, where σ^t2 is obtained recursively through (12). Afterwards, upon observing yl,m+1≤l≤m+n, the monitoring procedure is conducted by computing the CUSUM test statistic as in (8), and we declare it out of control when it exceeds the prescribed critical value under the nominal level α∈(0,1). Remark 1. In our proposed monitoring procedure, we use the theoretically obtained critical value c, as illustrated by its notable performance in Section 4. However, if m is not large enough relative to n, there is a chance that the monitoring procedure might be undermined by the parameter estimation. In this case, one may be able to obtain c empirically through a wild bootstrap approach with the following steps [21]: 1.  Estimate ω,α,β with ω^,α^,β^ from training sample y1,…,ym; 2.  Estimate σt2 recursively with σ^t2=ω^+α^yt−12+β^σ^t−12 and some initial values y0 and σ^0; 3.  Generate iid standard normal random variables ηtb, t=1,…,n, b=1,…,B, and construct a bootstrap sample ytb=σ^tηtb; 4.  Based on ytb, t=1,…,n, estimate ω,α,β with ω^*,α^*,β^*, and calculate the bootstrapped residuals ϵ^tb*=ytb/σ^tb* with σ^tb* obtained recursively by σ^tb*2=ω^*+α^*yt−1b2+β^*σ^t−1b*2; 5.  Based on these residuals, construct the monitoring process T^nmax,b(k), k=1,…,n, b=1,…,B, similarly to T^nmax in (8) with W^kb analogously defined to W^k; 6.  Finally, the critical value c is determined as the 100α% upper quantile of T^nmax,b=max1≤k≤nT^nmax,b(k) for b=1,…,B. The critical value c can be obtained via a bootstrap method similar to that for the GARCH(1,1) model in (7). In this case, we use ytb=σ˜tηt. Based on the obtained bootstrap sample ytb, applying the SVR method again, we can get the conditional volatility estimates σ^t*b2=g^*(yt−1b2,σ˜t−1b2), where σ˜tb2 are proxies obtained from ytb’s. Then, the bootstrapped residuals are obtained as ϵ^tb*=ytb/σ^t*b. The critical value is obtained from these similarly to Steps (5) and (6) addressed above. Instead of σ^t*b2, alternatively one might be able to consider using σ˜t*b2:=g^(yt−1b2,σ˜t−1b2) to obtain the residuals ϵ^tb*:=ytb/σ˜t*b. 4. Simulation Experiments We assess the proposed SVR-GARCH monitoring process to measure the performance when the time series is simulated from linear or nonlinear variants of GARCH models, such as GARCH(p,q), asymmetric GARCH(p,q) (AGARCH), GJR-GARCH(p,q), and Box–Cox transformed threshold GARCH(p,q) (BCTT-GARCH), specified as follows: GARCH(p,q):yt=σtϵt,σt2=ω+∑i=1pαiyt−i2+∑j=1qβjσt−j2,αi≥0,βj≥0;AGARCH(p,q):yt=σtϵt,σt2=ω+α(yt−1−b)2+βσt−12,α≥0,β≥0;GJR-GARCH(p,q):yt=σtϵt,σt2=ω+∑i=1pα1,iyt−i+2+α2,iyt−i−2+∑j=1qβjσt−j2,yt+=max(yt,0),yt−=−min(yt,0),α1,i≥0,α2,i≥0,βj≥0; BCTT-GARCH(p,q):yt=σtϵt,σt2=ω+∑i=1pα1,i(yt−i+2)δ+α2,i(yt−i−2)δ+∑j=1qβjσt−j21/δ,yt+=max(yt,0),yt−=−min(yt,0),α1,i≥0,α2,i≥0,βj≥0, where the orders p,q are fixed as 1 and the errors ϵt are iid. N(0,1) are random variables. In implementation, we generate a time series of length m from each of the above models and fit the SVR-GARCH model to it as presented in Section 3. In this procedure, we utilize the first 0.7m time series as a training set, and the rest as a validation set, where x denotes the largest integer not exceeding x. Subsequently, to create the circumstance, where the dataset for monitoring (testing sample) is observed sequentially, we initially design the underlying parametric model with a change located at point 1