Year: 2021

Abstract:
Sparse additive models have shown promising flexibility and interpretability in processing time series data. However, existing methods usually assume the time series data to be stationary and the innovation is sampled from a Gaussian distribution. Both assumptions are too stringent for heavy-tailed and non-stationary time series data that frequently arise in practice, such as finance and medical fields. To address these problems, we propose an adaptive sparse Huber additive model for robust forecasting in both non-Gaussian data and (non)stationary data. In theory, the generalization bounds of our estimator are established for both stationary and non-stationary time series data, which are independent of the widely used mixing conditions in learning theory of dependent observations. Moreover, the error bound for non-stationary time series contains a discrepancy measure for the shifts of the data distributions over time. Such a discrepancy measure can be estimated empirically and used as a penalty in our method. Experimental results on both synthetic and real-world benchmark datasets validate the effectiveness of the proposed method.

text: Sparse additive models have shown promising flexibility and interpretability in processing time series data. However, existing methods usually assume the time series data to be stationary and the innovation is sampled from a Gaussian distribution. Both assumptions are too stringent for heavy-tailed and non-stationary time series data that frequently arise in practice, such as finance and medical fields. To address these problems, we propose an adaptive sparse Huber additive model for robust forecasting in both non-Gaussian data and (non)stationary data. In theory, the generalization bounds of our estimator are established for both stationary and non-stationary time series data, which are independent of the widely used mixing conditions in learning theory of dependent observations. Moreover, the error bound for non-stationary time series contains a discrepancy measure for the shifts of the data distributions over time. Such a discrepancy measure can be estimated empirically and used as a penalty in our method. Experimental results on both synthetic and real-world benchmark datasets validate the effectiveness of the proposed method.

heading: 1 INTRODUCTION
text: Additive model has become one of the most powerful tools for time series analysis due to the exemplary monograph (Stone, 1985; Hastie & Tibshirani, 1990) and companion software (Chambers & Hastie, 1992). For the past two decades, the growing importance of algorithmic flexibility and interpretability motivates the development of various additive models along with theoretical explorations (Huang & Yang, 2004; Wang & Yang, 2007; Chu & Glymour, 2008; Song & Yang, 2010; Yang et al., 2018) and practical applications (Dominici et al., 2002; Wang & Brown, 2011; Ravindra et al., 2019; Bussmann et al., 2020). Although these aforementioned works have shown promising behaviours, the methods proposed in these works require some stringent assumptions on the stochastic process, e.g., various mixing conditions (Doukhan, 1994), stationary distribution and Gaussian innovation.
A number of attempts have been made to relax such stringent assumptions. For the purpose of dealing with non-Gaussian innovation, Qiu et al. (2015) develop an elliptical vector autoregressive model for estimating heavy-tailed stationary processes with parametric convergence analysis that reduces the influence of non-Gaussian innovation. Moreover, under stationarity and decaying β-mixing condition, Wong et al. (2020) derive nonasymptotic estimation error of lasso without assuming special parametric form of the data generating process. Stationarity and various mixing conditions are commonly adopted in many previous studies, see, e.g., (Mohri & Rostamizadeh, 2009; 2010; Kock & Callot, 2015; Wong et al., 2020). However, the mixing and stationary conditions are too stringent and not always valid in practice (Baillie, 1996; Kuznetsov & Mohri, 2020a). To relax the mixing and stationary conditions, Adams & Nobel (2010) prove asymptotic guarantees for stationary ergodic sequences. Agarwal & Duchi (2013) establish generalization bounds for asymptotically stationary (mixing) processes in the case of stable on-line learning algorithms. Kuznetsov & Mohri (2014) establish learning guarantees for fully non-stationary and mixing processes. Recently, Kuznetsov & Mohri (2020a) provide data-dependent generalization risk bounds for non-stationary and non-mixing stochastic processes.
However, both theoretically and practically, the exploration on robust additive models for nonstationary time series analysis is still limited. In this paper, we propose a class of sparse Huber additive models with theoretical guarantees. Our main contributions are summarized as follows:
• Algorithm design and theoretical guarantees: In Section 3.1, we first propose a novel sparse Huber additive model (SpHAM) by integrating Huber loss and sparsity-inducing `2,1-norm regularizer into an additive data-dependent hypothesis space. This proposed method can achieve robust forecasting and satisfactory inference (e.g., Granger causal discovery) simultaneously. In theory, Section 3.2 establishes the upper bound of the function approximation error of SpHAM by developing error decomposition technique (Wu et al., 2006; Chen et al., 2020) and employing sequential Rademacher complexity (Rakhlin et al., 2010; 2015; Kuznetsov & Mohri, 2020a). With properly selected scale parameters, the theoretical findings indicate that: a) For stationary time series, the function approximation consistency with convergence rate O(n− 1 2 ) can be pursued, even if the innovation is non-
Gaussian distribution (see Theorem 1 and Corollary 1 for more details). Moreover, this consistency analysis appears to be novel because no explicit data dependence assumptions are not imposed here, e.g., various mixing conditions used in (Doukhan, 1994; Mohri & Rostamizadeh, 2010; Zou et al., 2009; Wong et al., 2020); b) For non-stationary time series, the function approximation error is bounded by a discrepancy measure, which characterizes the drifts of the data distributions along with time (see Theorem 2 for more details). By penalizing such a discrepancy measure, Section 3.3 further proposes an adaptive SpHAM for dealing with non-stationary time series and provides its theoretical upper bound correspondingly (see Theorem 3 for more details).
• Optimization and empirical evaluations: The proposed SpHAM and adptive SpHAM can be implemented efficiently by Difference of Convex programming (DC programming) (Tao & An, 1998) and Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) (Beck & Teboulle, 2009). Experimental results on both synthetic and real-world benchmark CauseMe (Runge et al., 2019) validate the effectiveness of the proposed method.
Related works: There exist some interesting studies towards sparse additive model for time series analysis from both theoretical and practical viewpoints; see, e.g., (Chu & Glymour, 2008; Song & Yang, 2010; Yang et al., 2018). Although they show promising interpretability and modeling capacity, all of them require Gaussian innovation, stationarity and various mixing conditions. Moreover, these assumptions are not always valid in heavy-tailed and non-stationary time series data. In contrast to these previous additive models, we are seeking to formulate our method and investigate its asymptotic properties without resorting to such strict assumptions.
Robust forecasting is a vibrant area of research in time series analysis (Qiu et al., 2015; Wong et al., 2020). As one of the triumphs and milestones of robust statistics, the success stories of Huber algorithms are mainly on robust prediction tasks with i.i.d datasets (Huber, 1964; Huber & Ronchetti, 2009; Loh, 2017; Feng & Wu, 2020), and their extensions to time series prediction are fairly sparse. To our best knowledge, this is the first work that considers Huber additive models for non-stationary and dependent time series data.
Our theoretical analyses are inspired from the successful usage of sequential Rademacher complexity in Kuznetsov & Mohri (2020b;a). The generalization risk bounds are established in their works for a general scenario of non-stationary and non-mixing stochastic processes. However, it should be mentioned that we are chiefly concerned about the function approximation analysis, which is essentially different from theirs and is very crucial for Huber algorithms, since the convergence of generalization risk cannot imply the convergence of function approximation (Sun et al., 2019; Feng & Wu, 2020). Moreover, the development of analysis techniques (e.g., error decomposition and sequential Rademacher complexity) for assessing our method may shed light on other robust models for nonstationary time series analysis.
To better highlight the novelty of our method, Table 1 summarizes the algorithmic properties of our method and other related works, e.g., Sparse additive models for time series (TS-SpAM) (Yang et al., 2018), Discrepancy-based Forecasting (DBF) (Kuznetsov & Mohri, 2020a) and Robust Dantzigselector-type estimator (R-Dantzig) (Qiu et al., 2015).
The remainder of this paper is organized as follows. Section 2 recalls the background of additive models. Section 3 mainly provides our methods and theoretical guarantees. We provide empirical evaluations in Section 4. Finally, Section 5 concludes this paper.

heading: 2 PRELIMINARY
text: Let {Zt}∞t=−∞ be a stochastic time series with time index t, where variable Zt = (Xt, Y t) takes values in the compact input space X ⊂ Rp and the output space Y ⊂ R. We consider a common nonparametric model
Y t = f∗(Xt) + εt, E(εt) = 0, (1)
where f∗(·) is the ground truth function, and the innovation εt is i.i.d. across time t ∈ Z. For the sake of simplicity, we denote ρt and ρtX as the jointed distribution of (X
t, Y t) and the corresponding marginal distribution with respect to Xt, respectively. This setup actually covers a large number of scenarios commonly used in practice. For instance, the case Xt that contains p lagged values of Y t (e.g., Xt = Y (t−p)×· · ·×Y (t−1)) corresponds to the p-order autoregressive models. Moreover, this case can be viewed as a vector autoregressive model in the sense that input Xt includes the historical information of multiple variables.
Although such a nonparametric model (1) makes very few assumptions on data generation, the related nonparametric algorithms suffer so-called “curse of dimensionality”, see Fan & Gijbels (1996) for further discussion. An effective strategy for solving this problem is additive model. Usually, the additive structure is obtained by decomposing the input space X ∈ Rp into X = X1 × ... × Xp. Under the assumption that the ground truth admits an additive structure f∗ = ∑p j=1 f ∗ j , the additive model can be defined as
Y t = f∗1 (X 1 1 ) + · · ·+ f∗p (Xtp) + εt, (2)
where each component f∗j : Xj → R is a smooth function. In linear time series analysis, a weak stationarity condition (i.e., the first two moments of time series are time invariant) is preferred (Han et al., 2015; Qiu et al., 2015). In contrast, strict stationarity is primarily used if our focus is on nonlinear relationships (Fan & Yao, 2005). Definition 1. A stochastic process {Zt}∞t=−∞ is strictly stationary if (Z1, ..., Zt) and (Z1+k, ..., Zt+k) have the same joint distributions for any t ∈ Z and k ∈ Z.
Note that, if not otherwise stated, the stationarity in this paper refers to strict stationarity. Suppose that we are given T size time series data {(xt, yt)}Tt=1 ∈ ZT which are drawn from an additive data-generating model (2). Under stationarity condition and zero-mean Gaussian innovation with finite variance, the most existing methods that learn the ground truth f∗ usually integrate squared loss and a smoothness- or sparsity-inducing regularizer Ω(·) into a structural risk minimization scheme:
min f∈H T∑ t=1 (yt − p∑ j=1 fj(x t j)) 2 + Ω(f),
whereH := {f1 + ...+ fp : fj ∈ Hj , j = 1, ..., p} is an additive hypothesis space. Commonly, each subspaceHj can be reproducing kernel Hilbert space (Chen et al., 2017; Kandasamy & Yu, 2016; Raskutti et al., 2012), the orthogonal basis inducing space (Ravikumar et al., 2009; Meier et al., 2009;
Yang et al., 2018), or the composite function space with the neural network as a typical (Agarwal et al., 2020; Bussmann et al., 2020).
However, when facing heavy-tailed innovation, these methods may have degraded performance due to the amplification of the squared loss to large residuals. In statistic learning community, as one commonly used robust statistic, Huber loss is defined as
`σ(f(x t)− yt) = { (f(xt)− yt)2, if |f(xt)− yt| < σ 2σ|f(xt)− yt| − σ2, if |f(xt)− yt| ≥ σ, (3)
where σ is a positive hyper-parameter. Note that in the previous studies (Huber & Ronchetti, 2009; Loh, 2017), the hyper-parameter σ is set to be fixed according to the 95% asymptotic efficiency rule. However, Huber regression with a fixed scale parameter may not be able to learn the ground truth when the noise is asymmetric, as argued recently in Feng & Wu (2020); Sun et al. (2019). In this paper, we choose the scale parameter σ by relating it to the moment condition of the noise distribution and the sample size so that the resulting regression estimator can asymptotically converge to the ground truth function.

heading: 3 METHOD

heading: 3.1 SPARSE HUBER ADDITIVE MODELS
text: In this paper, we choose reproducing kernel Hilbert space (RKHS) HKj , j = 1, .., p, to form the the additive hypothesis space H, where each HKj is associated with a symmetric and positive semi-definite Mercer kernel Kj : Xj ×Xj → R. An additive RKHS is defined as
HK = {f1 + ...+ fp : fj ∈ HKj , j = 1, ..., p} (4) with kernel norm ‖f‖2K = inf{‖f‖2K1 + ...+ ‖f‖ 2 Kp }. By integrating the Huber loss (3), additive RKHS and kernel norm inducing regularizer into a Tikhonov regularization scheme, the kernel-based Huber additive model can be formulated as
f̂η = p∑ j=1 f̂η,j = arg min f= ∑p j=1 fj ,fj∈HKj { T∑ t=1 `σ(y t − p∑ j=1 fj(x t j)) + η p∑ j=1 τj‖fj‖2Kj}, (5)
where η is positive regularization parameter and τj is the weight for j-th kernel norm. The representer theorem (Wahba, 1990) ensures that f̂η can be represented as
f̂η = p∑ j=1 T∑ t=1 αηtjKj(x t j , ·), α η tj ∈ R.
To offer the method with sparsity, we consider the following sparsity-inducing penalty
Ω(f) := inf{ p∑ j=1 τj‖αj‖2 : f = p∑ j=1 T∑ t=1 αtjKj(x t j , ·)}.
LetHZ be an additive data dependent hypothesis space defined by
HZ = {f = p∑ j=1 T∑ t=1 αtjKj(x t j , ·) : αtj ∈ R}.
Then the SpHAM can be formulated as
f̂ = arg min f∈HZ { 1 T T∑ t=1 `σ(y t − p∑ j=1 fj(x t j)) + λΩ(f)}.
Denote Ktj = (Kj(x 1 j , x t j), ...,Kj(x T j , x t j)) ′ ∈ RT , α = (α′1, ..., α′p)′ ∈ RTp and αj = (α1j , ..., αTj) ′ ∈ RT , where α′j here refers to the transposition of αj for avoiding confusion with time T . The SpHAM can be represented as
f̂ = p∑ j=1 T∑ t=1 αλtjKj(x t j , ·) (6)
with
αλ = arg min αj∈RT ,j=1,...,p { 1 T T∑ t=1 `σ(y t − p∑ j=1 (Ktj) ′αj) + λ p∑ j=1 τj‖αj‖2}. (7)
The optimization problem (7) can be solved by Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) (Beck & Teboulle, 2009). We provide the detailed optimization procedure in Appendix E due to space limitation.
Remark 1. Inspired by the studies on group additive models (Yin et al., 2012; Pan & Zhu, 2017; Chen et al., 2017), the modeling capacity of our proposed method can be improved by embedding the (Hierarchical) group structure. For instance, group SpHAM can be easily formulated by replacing the direct decomposition {Xj}pj=1,Xj ∈ R with the subgroups decomposition {Xd}Dd=1, where each Xd ⊂ X may cover more than one variable.

heading: 3.2 ASYMPTOTIC THEORY ANALYSIS
text: In this section, we provide the theoretical analysis of SpHAM. It is natural and necessary to investigate the following three problems:
• The function estimation consistency and convergence rate under stationary time series data (See Theorem 1 in Section 3.2.1);
• The upper bound of the function estimation error under non-stationary time series data (See Theorem 2 in Section 3.2.2);
• Whether the theoretical findings can inspire an adaptive SpHAM for dealing with nonstationary time series (See Theorem 3 in Section 3.3)
Due to space limitation, the high-level outline and detailed proofs are provided in the Appendix A-D.

heading: 3.2.1 FUNCTION APPROXIMATION ANALYSIS FOR STATIONARY TIME SERIES
text: Through this paper, the marginal distribution w.r.t Y t is assumed to be almost everywhere supported on [−M,M ] for some M ≥ 0. Assumption 1. Assume that |Y t|, ∀t ∈ Z, is bounded and there exists a constant c > 0 such that E|Y t|1+c <∞, ∀t ∈ Z.
The moment condition in Assumption 1 is rather weak in the sense that the response variable Y t possesses infinite variance. The same condition also applies to the distributions of the innovation εt, implying that heavy-tailed innovation is allowed.
Assumption 2. Let κ = supx∈X √ Kj(x, x) <∞,∀j = 1, ..., p.
Assumption 2 only requires that the kernel is bounded under compact X , which holds for all Mercer kernels (e.g., Gaussian kernel) and has been used in many literatures; see, e.g., (Wu et al., 2006; Steinwart & Christmann, 2008; Wu & Zhou, 2008).
The following definitions are needed for Assumption 3. For any j = 1, ..., p, we define a kernel integral operator LKj ,T+1 : L2(ρ T+1 Xj )→ L2(ρ T+1 Xj ) associated with the kernel Kj by
LKj ,T+1(f)(x T+1 j ) = ∫ Xj Kj(x T+1 j , uj)f(uj)dρ T+1 Xj (uj).
Note thatLKj ,T+1 is a compact and positive operator onL2(ρ T+1 Xj ). According to Mercer theorem, we can find the corresponding normalized eigenpairs {(ζji , ψ j i )}i≥1 such that {ψ j i }i≥1 is an orthonormal basis of L2(ρT+1X ) and ζ j i → 0 as i→∞. Then for given r > 0, we defined the r-th power LrKj ,T+1 by LrKj ,T+1(
∑ i≥1 βjiψ j i ) = ∑ i≥1 βji (ζ j i ) rψji .
Assumption 3. We assume that f∗j : Xj → R,∀j = 1, ..., p is a function of the form f∗j = LrKj ,T+1(g ∗ j ),∀r ∈ (0, 12 ] with some g ∗ j ∈ L2(ρ T+1 Xj ),∀h ∈ Z.
Assumption 3 is a natural extension from i.i.d setting (Assumption 1 in (Wu et al., 2006; Christmann & Zhou, 2016)) to non-i.i.d time series. Indeed, this assumption stands in most practical cases. For instance, if r = 0.5, the ground truth function f∗ needs to be a real-valued function in the RKHS HK . The RKHS admits a large class of bounded and continuous real-valued functions that are in Hilbert space. Moreover, this assumption has been widely used in learning theory; please refer to (Smale & Zhou, 2003; Cucker & Zhou, 2007) for more practical discussions (e.g., 0 < r < 0.5). Theorem 1. Suppose that the process {Zt}∞t=−∞ is stationary. Let Assumptions 1-3 be true. By taking σ = Tm, η = T− 1 4r and λ = T− 1 4r−m, we have for any 0 < δ < 1
‖f̂ − f∗‖2 L2(ρ T+1 X ) ≤ C̃ log(1/δ)TΨ(m,c,r)
with confidence at least 1− δ, where C̃ is a positive constant independently of T, λ, η, δ and σ and
Ψ(m, c, r) =
{ max{− 12 ,m− 1,−cm+ 1 4r +m− 1}, if m ≤ 1− 1 4r
max{− 12 ,m− 1,−cm+ 1 2r + 2m− 2}, if m > 1− 1 4r .
Remark 2. In the stationary case, our bound appears to be novel for the following reasons: a) the result is completely independent of the various mixing conditions which are widely used in the theory analysis of non-i.i.d dependent time series (Mohri & Rostamizadeh, 2009; Yang et al., 2018; Mohri & Rostamizadeh, 2010; Guo & Shi, 2011; Qiu et al., 2015); b) compared with Yang et al. (2018), the innovation assumption is rather weak in the sense that the innovation possesses infinite variance and thus admits a heavy-tailed distribution. Corollary 1. Suppose that the process {Zt}∞t=−∞ is stationary. Let all the conditions in Theorem 1 be true. We then have for any 0 < δ < 1
‖f̂ − f∗‖2 L2(ρ T+1 X ) ≤ C̃ log(1/δ)TΨ(m,c)
with confidence at least 1− δ, where
Ψ(m, c) =
{ max{− 12 ,−cm+m− 1 2}, if m ≤ 1 2
max{m− 1,−cm+ +2m− 1}, if m > 12 .
Figure 1 summaries the convergence rates in Corollary 1 by taking different σ and c. If the innovation is Gaussian distribution with finite variance (i.e., Assumption 1 holds for any c > 1), one can arbitrarily select a σ = Tm (0 < m < 12 ) to obtain convergence rates O(n− 1 2 ). Moreover, we can see that the convergence rate will decrease as m increases. Combined with the conclusion in Lemma 1 (i.e., the equivalence relation between Huber loss based empirical risk and MSE as σ → ∞), it indicates that σ indeed plays a tradeoff role between algorithmic robustness and variance-reduction. For the weak moment condition (e.g., 0 < c < 1), one may get slower convergence rates (e.g., O(n(1−c)m− 1 2 ) or
O(n(2−c)m−1)), which also coincides with our intuitive understanding that small σ may be conducive to robust forecasting. Note that our method will not converge when σ and c are both located in the white area in Figure 1.
As a comparison, due to the non-robustness of the squared loss, the most existing convergence rates are established under the assumption that the innovation is Gaussian distribution with finite variance (see, e.g., Yang et al. (2018); Han et al. (2015); Wang & Yang (2007); Kock & Callot (2015)). However, from Theorem 1 and Corollary 1, the asymptotic convergence of SpHAM can be obtained under weaker moment condition, which verifies the robustness of our method. Recall the learning rate O(n 2d 2d+1 ) derived in Yang et al. (2018), where d is the order of smoothness of the component function fj , j = 1, ..., p. Relatively slow convergence rate O(n− 1 2 ) we obtained indicates the sacrifice for the absence of mixing condition.

heading: 3.2.2 FUNCTION APPROXIMATION ANALYSIS FOR NON-STATIONARY TIME SERIES
text: In nonstationary time series setting, different Zts may follow different distributions. Thus we introduce the weights s = {s1, ..., sT } and assign them to the losses made on different sample points in terms of their relevance to forecasting the future ZT+1. For given weights {st}Tt=1, we define the estimator f̂ s = ∑p j=1 f̂ s j as the minimizer of the following weighted objective
Esλ(f̂ s) := min f∈HZ { T∑ t=1 st`σ(y t − p∑ j=1 fj(x t j)) + λΩ(f)}. (8)
Next, we introduce a discrepancy measure to characterize the discrepancy between the target distribution and the distributions of observations (Kuznetsov & Mohri, 2020a).
Definition 2. For any f ∈ HZ , the discrepancy measure with respect to Huber loss is defined as
disc(s) := sup f∈HZ
{ E`σ(f(xT+1)− yT+1)− T∑ t=1 stE`σ(f(xt)− yt) } .
The discrepancy measure is a natural measure of the non-stationarity of the stochastic process {Zt}∞t=−∞ with respect to both the loss function `σ and the hypothesis setHZ . Theorem 2. Let Assumptions 1-3 be true. We assume that f∗j ∈ HKj ,∀j = 1, ..., p. By taking σ = T 1 2c , λ = T−1 and η = T− 1 2 , we have for any 0 < δ < 1
‖f̂ s − f∗‖2 L2(ρ T+1 X ) ≤ disc(s) + Esλ(f̂ s) + C̃1‖s‖2T 1 2c + C̃2 log(1/δ)T − 12
with confidence at least 1− δ, where C̃1, C̃2 are two positive constants independently of T, λ, η, δ and σ.
Note that there are several existing studies towards analyzing nonstationary and non-mixing time series, see, e.g., Kuznetsov & Mohri (2020a). Different from them, the error bound we derived is with respect to the function approximation rather than generalization risk, which is very crucial for Huber regression problem, since the convergence of generalization risk cannot imply the convergence of function approximation (Sun et al., 2019; Feng & Wu, 2020).

heading: 3.3 ADAPTIVE SPARSE HUBER ADDITIVE MODEL FOR NONSTATIONARY TIME SERIES
text: Theorem 2 illustrates that we shall minimize the following optimization problem for non-stationary time series forecasting:
min f∈HZ { T∑ t=1 st`σ(yt − p∑ j=1 fj(x t j)) + disc(s) + λ1Ω(f) + λ2T 1 2c ‖s‖2},
where λ1 and λ2 are two positive regularization parameters, and c is a positive constant introduced in Assumption 1. Although the discrepancy measure disc(s) is crucial for such an optimization problem, we cannot obtain its exact value since we do not have access to the distributions of Zt, t ∈ Z. Hence, we need to estimate the approximated discrepancy from given data. Inspired by (Kuznetsov & Mohri, 2020a), one natural and necessary assumption is that there is an underlying representation relationship between distribution ρT+1 and distributions ρt, t = 1, ..., T .
Assumption 4. Denote by a probability set q∗ = {q∗t }Tt=1 with ∑T t=1 q ∗ t = 1. We assume that the following term is sufficiently small:
disc(q∗) := sup f∈HZ [E`σ(f(xT+1)− yT+1)− T∑ t=1 qtE`σ(f(xt)− yt)].
Note that the priori q∗ can be any distribution, which shall be given empirically according to the trend of the time series. For instance, in a particular scenario of the distribution ρT+1 does not change drastically compared with the distributions ρt, t = 1, ..., T , Kuznetsov & Mohri (2020a) have proven that the Assumption 4 holds if the s∗ is an uniform distribution over last l > 0 observations, where l is a hyper-parameter that can be tuned in practical applications. Theorem 3. Let Assumption 4 and the conditions in Theorem 2 be true. We have for any 0 < δ < 1
‖f̂s − f∗‖2 L2(ρ T+1 X ) ≤ sup f∈HZ T∑ t=1 (q∗t − st)`σ(f(xt)− yt) + Esλ(f̂ s)
+C̃1 log(1/δ)(‖q∗ − s‖2 + ‖s‖2T 1 2c ) + C̃2 log(1/δ)T − 12
with confidence at least 1− δ, where C̃1, C̃2 are two constants independently of T, λ, η, δ and σ.
Finally, the optimization problem of adaptive SpHAM can be formulated as following two stages:

heading: Step A: finding the weight ŝ:
text: ŝ = arg min s { sup f∈HZ T∑ t=1 (q∗t − st)`σ(f(xt)− yt) + λ1‖q∗ − s‖22 + λ2T 1 2c ‖s‖22} (9)
Step B: forecasting:
f̂ s(xT+1) = p∑ j=1 T∑ t=1 αstjK(x t j , x T+1 j ),
where
αs = arg min αj∈RT ,j=1,...,p { T∑ t=1 ŝt`σ(y t − p∑ j=1 (Ktj) ′αj) + λ p∑ j=1 τj‖αj‖2}. (10)
The optimization problems (9) can be solved by standard gradient descent, where, for each step, we employ DC programming (Tao & An, 1998) to evaluate the supremum over f ∈ HZ . Similarly to the strategy for optimization problem (7), we use Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) (Beck & Teboulle, 2009) for Step B. We provide the detailed procedure in Appendix E.

heading: 4 EXPERIMENT
text: This section validates the effectiveness of SpHAM and adaptive SpHAM. In all experiments, the Gaussian kernel Kj(u, v) = exp(−‖u−v‖ 2 2
2d2 ), where j = 1, ..., p and bandwidth d > 0, is employed for constructing the additive data dependent hypothesis space. Due to limited space, we only present the evaluations of part of the synthetic time series data. The real-world experiments can be found in Appendix F. In the following, we consider three synthetic examples as below:
Example A: Inspired by (Kuznetsov & Mohri, 2020a), a stationary time series is generated according to the non-linear additive autoregressive model:
Y t = 3
2 sin(
π 2 Y t−2)− sin(π 2 Y t−3) + εt,
where the innovation εts are i.i.d. drawn from Gaussian distributionN(0, 0.1) and Student distribution with freedom 2, respectively.
Example B: A stochastic time series with interaction term is generated by
Y t = sin(Y t−1Y t−2) + t.
Example B considers the data generating process with interaction term Y t−1Y t−2. In order to improve the modeling capacity of our SpHAM, according to Remark 1, Group SpHAM is introduced here by dividing the input space into 3 subgroups {X1}, {X2} and {X1,X2}. Example C: Inspired by (Kuznetsov & Mohri, 2020a), a time series with smooth drift is generated by
Y t = t
400 sin(Y t−1) + t.
Hyper-parameter selection and evaluation criterions: Recall that SpHAM algorithm requires three hyper-parameters: regularization parameter λ, bandwidth of kernel d and Huber parameter σ. For Huber parameter σ selection, we set σ = T 1 48 based on the suggestion in Theorem 1 and Corollary 1. Moreover, we grid search the hyper-parameters: the bandwidth d ∈ {0.1, 0.5, 1, 5}, the regularization parameters λ, λ1, λ2 ∈ {0.001, 0.01, 0.1, 1, 10, 100, 1000} and l ∈ {100, 150, 200, ..., 350}. The evaluation criterions for forecasting used here contains Average Sample Error(ASE)= 1N √∑N t=1(f̂(x t)− yt)2 and True Deviation (TD)= 1N √∑N t=1(f̂(x t)− f∗(xt))2, where N is the number of test samples.
For each example, we generate time series with 4000 sample points. The samples at time t = {1500, 1501, ..., 1899} are used as a training set, and the samples at time {1900, 1901, 1902} are considered as the test set. The competitors include standard Autoregression (AR), TS-SpAM (Yang et al., 2018), Gaussian processes regression (GP) and Vanilla Long Short-Term Memory (LSTM). The tuning parameters of other competing methods such as LSTM and GP, are chosen according to their original sklearn package in python (Pedregosa et al., 2011). All the evaluations are repeated for 50 times. The average results for Examples A-C are presented in Table 2. Moreover, in Appendix F.1, we investigate the impact of the σ on the performance of our method.
Table 2 verifies that robust forecasting of our SpHAM in presence of heavy-tailed t noise. Moreover, when facing the stochastic time series with interaction term, the result (Example B) indicates the modelling capacity of our SpHAM can be further improved by embedding group structure. The result (Example C) shows the promising performance of additive SpHAM on nonstationary time series forecasting.

heading: 5 CONCLUSION
text: We propose an adaptive sparse Huber additive model by integrating Huber loss and `2,1-norm regularizer into an additive data dependent hypothesis space. We theoretically explore the asymptotic properties of our method for both non-Gaussian and (non)stationary time series. Experimental results on both synthetic and real-world data validate the effectiveness of the proposed method.

heading: A PROOF SKETCH
text: This section sketches the proof focusing on the conceptual aspects. The full proofs are provided in Appendix B-D.
Figure 2 summaries the important ingredients for the proof of Theorem 1. Lemma 1 describes the properties off̂ and Lemma 2 illustrates the relation between huber-based risk and MSE. Following these two lemmas, we then decompose the function error into three important parts: sample error, hypothesis error and approximation error. To bound the sample error, we need to employ the sequential Rademacher complexity (Kuznetsov & Mohri, 2020a), which is the key to measuring the capacity of the data-dependent hypothesis space for non-i.i.d data. Lemmas 3-4 give an important concentration inequality and the upper bound of sequential Rademacher complexity, respectively. Finally, Theorem 1 is obtained by combining Propositions 2-4. The proof of Theorems 2 is similar to the above process, except that we develop a different error decomposition in Proposition 5. On the basis of Theorem 2, by applying an assumption on data distribution, we derive Theorem 3.
The full proofs of Theorems 1-3 are provided in Sections C-E, respectively.

heading: B PROOF OF THEOREM 1
text: We first illustrate two key properties of f̂ .
Lemma 1. Let Assumptions 1-2 be true. From the definition of f̂ in Eq. (6), there hold
‖f̂‖∞ ≤Mf̂ := 2κ2M2
λT minj=1,...,p τj .
and
Ω(f̂) = p∑ j=1 τj‖αλj ‖2 ≤ κ−2Mf̂ = 2M2 λT minj=1,...,p τj ,
where M is a positive constant such that |Y t| ≤M,∀t ∈ Z. Proof. Denote by EσT (f) := ∑T t=1 `σ(f(x
(t))−y(t)) for notational convenience. From the definition of f̂ , we know that EσT (f̂) + λΩ(f̂) ≤ EσT (0). It means that
Ω(f̂) = p∑ j=1 τj‖αλj ‖2 ≤ λ−1EσT (0).
From the definition of Huber loss (3), we obtain
EσT (0) = 1
T T∑ t=1 `σ(y (t)) ≤ σ2, if |y(t)| ≤ σ, ∀t = 1, ..., T,
and
EσT (0) = 1
T T∑ t=1 `σ(y (t)) = 2σ|y(t)| − σ2 ≤ 2σM − σ2, if |y(t)| > σ, ∀t = 1, ..., T.
Therefore, p∑ j=1 τj‖αλj ‖2 ≤ 2M2 λT minj=1,...,p τj .
According to the property of RKHS, we conclude that
‖f̂‖∞ ≤ κ‖f̂‖K ≤ κ √√√√ p∑ j=1 ‖f̂j‖2Kj ≤ κ 2 p∑ j=1 ‖αλj ‖2 ≤ 2κ2M2 λT minj=1,...,p τj .
This completes the proof by denoting Mf̂ = 2κ2M2 λT minj=1,...,p τj .
Let H̄Z be a data-dependent hypothesis space constrained by sparsity-inducing regularizer Ω(·), i.e.,
H̄Z = {f(x) = p∑ j=1 T∑ i=1 αijKj(x (t) j , xj) : p∑ j=1 τj‖αj‖2 ≤ 2M2 λT minj=1,..,p τj }.
For any h ∈ Z, denote by EσT+1(f) := E`σ(f(xT+1)−yT+1) and ET+1(f) := E(f(xT+1)−yT+1)2, respectively. In the following, we describe the relationship between EσT+1(f) and ET+1(f), ∀f ∈ H̄Z . Note that a similar proof is given in (Feng & Wu, 2020) for the i.i.d case. Lemma 2. Let Assumptions 1-2 be true. For any h ∈ Z and f ∈ H̄Z , there holds∣∣∣[EσT+1(f)− EσT+1(f∗)]− [ET+1(f)− ET+1(f∗)]∣∣∣ ≤ Mcσc , where both c and Mc = 23+c(Mf∗ +Mf̂ + 1) 2E(|Y T+1|1+c) are positive constants.
Proof. For any σ > max{M +Mf̂ , 1}, we denote two events IY and IIY as follows
IY T+1 = {Y T+1 : |Y T+1| ≥ σ/2} and IIY T+1 = {Y T+1 : |Y T+1| < σ/2}. Then, for any f ∈ H̄Z and ‖f∗‖∞ = Mf∗ <∞, we have∣∣∣[EσT+1(f)− EσT+1(f∗)]− [ET+1(f)− ET+1(f∗)]∣∣∣
= ∣∣∣[ ∫ Z `σ(f(x T+1)− yT+1)dρT+1 − ∫ Z `σ(f ∗(xT+1)− yT+1)dρT+1 ] − ‖f − f∗‖2 L2(ρ T+1 X ) ∣∣∣ =
∣∣∣[ ∫ X ∫ Y `σ(f(x T+1)− yT+1)− `σ(f∗(xT+1)− yT+1)dρT+1Y|X dρ T+1 X ] − ‖f − f∗‖2 L2(ρ T+1 X ) ∣∣∣ ≤
∣∣∣ ∫ X ∫ I Y T+1 ∪II Y T+1 `σ(f(x T+1)− yT+1)− `σ(f∗(xT+1)− yT+1)dρT+1Y|X dρ T+1 X ∣∣∣ + ∣∣∣ ∫ X ∫ I Y T+1 ∪II Y T+1 (f(xT+1)− yT+1)2 − (f∗(xT+1)− yT+1)2dρT+1Y|X dρ T+1 X
∣∣∣ ≤
∣∣∣ ∫ X ∫ I Y T+1 `σ(f(x T+1)− yT+1)− `σ(f∗(xT+1)− yT+1)dρT+1Y|X dρ T+1 X ∣∣∣ + ∣∣∣ ∫ X ∫ I Y T+1 (f(xT+1)− yT+1)2 − (f∗(xT+1)− yT+1)2dρT+1Y|X dρ T+1 X ∣∣∣.
The last inequality is based on the fact that
|yT+1 − f(xT+1)| < |yT+1|+ ‖f‖∞ < σ, ∀(xT+1, yT+1) ∈ X × IIY T+1 and hence
`σ(f(x T+1)− yT+1) = (f(xT+1)− yT+1)2.
Recalling that the Huber loss is 2σ-Lipschitz continuous, then we get∣∣∣ ∫ X ∫ IY T+1 `σ(f(x T+1)− yT+1)− `σ(f∗(xT+1)− yT+1)dρT+1Y|X dρ T+1 X ∣∣∣ ≤ 2σ
∣∣∣ ∫ X ∫ IY T+1 |f(x)− f∗(x)|dρT+1Y|X dρ T+1 X ∣∣∣ ≤ 2σ‖f − f∗‖∞Prob{IY }.
According to Markov’s inequality, for any constant c > 0, we further have
Prob{IY T+1} = Prob{|Y T+1| ≥ σ/2} = Prob{|Y T+1|1+c ≥ (σ/2)1+c} ≤ 21+cE(|Y T+1|1+c)
σ1+c .
Then we have ∣∣∣ ∫ X ∫ IY T+1 `σ(f(x T+1)− yT+1)− `σ(f∗(xT+1)− yT+1)dρT+1Y|X dρ T+1 X ∣∣∣ ≤ 2
2+c‖f − f∗‖∞E(|Y T+1|1+c) σc .
Moreover, the second term in the right-hand side can be bounded by∣∣∣ ∫ X ∫ IY T+1 (f(xT+1)− yT+1)2 − (f∗(xT+1)− yT+1)2dρT+1Y|X dρ T+1 X ∣∣∣ ≤ ‖f − f∗‖∞ ∫ X ∫ IY T+1 |2yT+1 − f(xT+1)− f∗(xT+1)|dρT+1Y|X dρ T+1 X
≤ ‖f − f∗‖∞ [ ∫
IY T+1
2|yT+1|dρT+1Y + (‖f‖∞ + ‖f ∗‖∞)Prob(IY T+1)
] .
According to Holder inequality, we have∫ IY T+1 |Y T+1|dρT+1Y ≤ (Prob(IY T+1)) c 1+c (E|Y T+1|1+c) 1 1+c ≤ 2 cE(|Y T+1|1+c) σc .
Then we can deduce that ‖f − f∗‖∞ [ ∫
IY T+1
2|yT+1|dρT+1Y + (‖f‖∞ + ‖f ∗‖∞)Prob(IY T+1) ] ≤ (2
2+c + 21+c)‖f − f∗‖∞E(|Y T+1|1+c) σc + 21+c(‖f‖∞ + ‖f∗‖∞)2E(|Y T+1|1+c) σ1+c
Finally, there holds∣∣∣[ ∫ Z `σ(f(x T+1)− Y T+1)− ∫ Z `σ(f ∗(xT+1)− Y T+1) ] − ‖f − f∗‖2 L2(ρ T+1 X ) ∣∣∣ . Mc σc
with Mc = 2 3+c(Mf∗ +Mf̂ + 1) 2E(|Y T+1|1+c).
This completes the proof.
Define a stepping-stone function with respect to the distribution ρT+1:
f∗η,T+1 = p∑ j=1 f∗η,T+1,j = arg min f∈HK E`σ(f(xT+1)− Y T+1) + η p∑ j=1 τj‖fj‖2Kj .
To establish the bound of function approximation for stationary time series setting, we make the following error decomposition.
Proposition 1. Let Assumptions 1-2 be true. For any f ∈ H̄Z , there holds
‖f − f∗‖2 L2(ρ T+1 X ) ≤ E1 + E2 + E3 + 2Mc σc ,
where E1 = {EσT+1(f)− EσT (f)}+ {EσT (f∗η,T+1)− EσT+1(f∗η,T+1)}
E2 = ET+1(f∗η,T+1)− ET+1(f∗) + η p∑ j=1 τj‖f∗η,T+1,j‖2Kj
and
E3 = {EσT (f̂) + λΩ(f̂)− EσT (f̂η)− η p∑ j=1 τj‖f̂η,j‖2Kj}.
Proof. According to Lemma 2, for any f ∈ H̄Z , we can make following error decomposition
‖f − f∗‖2 L2(ρ T+1 X ) = ET+1(f)− ET+1(f∗) ≤ EσT+1(f)− EσT+1(f∗) + Mc σc
= {EσT+1(f)− EσT (f)}+ {EσT (f∗η,T+1)− EσT+1(f∗η,T+1)}︸ ︷︷ ︸ E1 +{EσT (f)− EσT (f∗η,T+1)}
+ EσT+1(f∗η,T+1)− EσT+1(f∗) + Mc σc
≤ E1 + ET+1(f∗η,T+1)− ET+1(f∗) + η p∑ j=1
τj‖f∗η,T+1,j‖2Kj︸ ︷︷ ︸ E2 +{EσT (f) + λΩ(f)− EσT (f∗η,T+1)
− η p∑ j=1 τj‖f∗η,T+1,j‖2Kj}+ 2Mc σc
≤ E1 + E2 + {EσT (f) + λΩ(f)− EσT (f̂η)− η p∑ j=1
τj‖f̂η,j‖2Kj}︸ ︷︷ ︸ E3 + 2Mc σc
≤ E1 + E2 + E3 + 2Mc σc .
In statistical machine learning community, we call E1, E2 and E3 sample error, approximation error and hypothesis error, respectively. The sample error E1 describes the divergence between the empirical risk EσT (f) and the expected risk Eσ(f). The hypothesis error E2 characterizes the difference between the empirical regularized risks with HK and HZ . The approximation error measures the approximation ability of RKHSHK toH.

heading: B.1 SAMPLE ERROR
text: In this section, we focus on providing the bound of sample error E1. Unfortunately, the traditional tools for complexity analysis such as covering number (Ron, 2000; Guo & Shi, 2011), Rademacher complexity (Mohri & Rostamizadeh, 2009)) and concentration inequalities ((Wu et al., 2007)), cannot be applied into this non-i.i.d. setting directly. To solve this problem, we employ the sequential Rademacher complexity developed in (Rakhlin et al., 2010; Kuznetsov & Mohri, 2020a). We define a function-based random variable as
ξf (z) = `σ(f(x)− y)− `σ(f∗η,T+1(x)− y), f ∈ H̄Z .
Now, we turn to establish the bound of
Eξf (zT+1)− 1
T T∑ t=1 ξf (z t), ∀f ∈ H̄Z .
A necessary ingredient needed for our analysis is data-dependent sequential Rademacher complexity Rakhlin et al. (2010), which we review in the following. We adopt the following definition of a complete binary tree: a Z-valued complete binary tree v is a sequence (v1, ..., vT ) of T mappings vt : {±1}t−1 → Z, t ∈ [1, T ]. A path in the tree is γ = (γ1, ..., γT−1) ∈ {±1}T−1. To simplify the notation, we will write vt(γ) instead vt(γ1, ..., γt−1), even though vt depends only on the first t− 1 elements of γ. The following definition generalizes the classical notion of Rademacher complexity to sequential setting.
Definition 3. The sequential Rademacher complexityRT (G) of a function class G is defined by
RT (G) = sup v E[sup g∈G
1
T T∑ t=1 γtg(vt(γ))],
where the supremum is taken over all complete binary trees of depth T with values in Z and γ is a sequence of Rademacher random variables.
Based on the definition of sequential Rademacher complexity, we have the following concentration inequality:
Lemma 3. Suppose that the time series {Zt}∞t=−∞ is strictly stationary. For any δ > 0, with probability at least 1− δ, the following inequality holds for all f ∈ H̄Z and all α > 0
Eξf (zT+1) ≤ 1
T T∑ t=1 ξf (z t) + 1√ T + 6M`σ √ 4π log TRT (G) + M`σ √ 8 log 1δ√ T ,
where G := {`σ(f, z) + `σ(fη, z) : f ∈ H̄Z} and M`σ is a positive constant such that `σ(f(x(t))− y(t)) ≤M`σ for any f ∈ H̄Z and t ∈ Z.
According to the property of sequential Rademacher complexity (see Proposition 14 in (Rakhlin et al., 2015)), we have
RT (G) = RT (`σ ◦ H̄Z , z) ≤ 16σ(1 + 4 √ 2 log3/2(eT 2))RT (H̄Z) (11)
with the 2σ-Lipschitz constant of `σ(·).
Lemma 4. Under Assumptions 1-2, there holds
RT (G) ≤ 32σ[1 + 4
√ 2 log3/2(eT 2)]M2p 1 2κ
T .
Proof. In fact, to every kernel Kj , we can associate a feature map φ with inner product < φ(xj), φ(x ′ j) >= K(xj , x ′ j) for any xj , x ′ j ∈ Xj and j = 1, ..., p. Then the constrained datadependent hypothesis space H̄Z can be rewritten as
H̄Z = {f(x) = p∑ j=1 T∑ t=1 αijφ(x (t) j )φ(xj) : p∑ j=1 τj‖αj‖2 ≤ 2M2 λT minj=1,..,p τj }.
By direct computation, we have
RT (H̄Z) = sup {x̄t,ȳt}Tt=1 Eσ[ sup α∈RTp
1
T T∑ t=1 σt( p∑ j=1 T∑ i=1 αijφ(x (i) j )φ(x̄ (t) j (σ))− y (t))]
= 1
T sup {x̄(t)}Tt=1 Eσ sup α∈RTp T∑ t=1 σt p∑ j=1 T∑ i=1 αijφ(x (i) j )φ(x̄ (t) j (σ))− 1 T sup {ȳ(t)}Tt=1 Eσ T∑ t=1 σty (t)(σ)
= 1
T sup {x̄(t)}Tt=1 Eσ sup α∈RTp T∑ t=1 σt p∑ j=1 T∑ i=1 αijφ(x (i) j )φ(x̄ (t) j (σ))
= 1
T sup {x̄(t)}Tt=1 Eσ sup α∈RTp T∑ t=1 σt p∑ j=1 T∑ i=1 αijK(x (i) j , x̄ (t) j (σ))
= 1
T sup {x̄(t)}Tt=1 Eσ sup α∈RTp α′ T∑ t=1 σtKt(σ), Kt(σ) = (K(x (1) 1 , x̄ (t) 1 (σ)), ...,K(x (T ) p , x̄ (t) p (σ)) ′ ∈ RTp
≤ 1 T sup α∈RTp ‖α′‖ sup {x̄t}Tt=1 Eσ‖ T∑ t=1 σtKt(σ)‖2
≤ 2M 2
T 2 sup {x̄(t)}Tt=1 Eσ‖ T∑ t=1 σtKt(σ)‖2
= 2M2
T 2 sup
{x̄(t)}Tt=1 √√√√Eσ[ T∑ s,t=1 σtσsK′t(σ)Ks(σ)]
≤ 2M 2
T 2 sup
{x̄(t)}Tt=1 √√√√Eσ[ T∑ t=1 K′t(σ)Kt(σ)]
≤ 2M 2p 1 2κ
T
Thus, combining above result with Eq. (11). We get the desirable result.
Combining the results in Lemma 4 with the inequality in Lemma 3, we then obtain the bound of sample error E1 for the strictly stationary time series setting. Proposition 2. Under Assumptions 1-2, for any δ > 0, the following inequality holds for all f ∈ H̄Z and all α > 0
E1 ≤ 1√ T
+ 6M`σ √ 4π log TR+ M`σ √ 8 log 1δ√ T
with probability at least 1− δ, where
R = 32σ[1 + 4
√ 2 log3/2(eT 2)]M2p 1 2κ
T .

heading: B.2 APPROXIMATION ERROR
text: Since the output function f∗η,T+1 is the optimal estimator in RKHSHK for time T + 1, the learning rates of the learning algorithm indeed depend on the approximation ability of the hypothesis spaceHK with respect to the optimal risk E(f∗) measured by the approximation error E2. For any j = 1, ..., p and h ∈ Z, we first define the kernel integral operator LKj ,T+1 : L2(ρT+1Xj )→ L2(ρ T+1 Xj ) associated with the kernel Kj by
LKj ,T+1(f)(x T+1 j ) = ∫ Xj Kj(x T+1 j , uj)f(uj)dρ T+1 Xj (uj),∀h ∈ Z.
Note that LKj ,T+1 is a compact and positive operator on L2(ρ T+1 Xj ). According to Mercer theorem, we can find the corresponding normalized eigenpairs {(ζjh,i, ψ j h,i)}i≥1 such that {ψ j h,i}i≥1 is an orthonormal basis of L2(ρT+1X ) and ζ j h,i → 0 as i→∞. Then for given r > 0, we defined the r-th power LrKj ,T+1 by
LrKj ,T+1( ∑ i≥1 βjh,iψ j h,i) = ∑ i≥1 βjh,i(ζ j h,i) rψjh,i.
We introduce an intermediate function as follows: f̃η,T+1,j = (LKj ,T+1 + ητjI) −1LKj ,T+1f ∗ j ,∀j = 1, ..., p.
Lemma 5. Under Assumption 3, for the intermediate function f̃η,T+1,j defined above, there holds
E‖f̃η,T+1,j(xT+1j )− f ∗ j (x T+1 j )‖+ ητj‖f̃η,T+1,j(x T+1 j )‖ 2 Kj ≤ 2(ητj) 2r‖L−rKj ,T+1f ∗ j ‖22.
Proof. Under Assumption 3, for any h ∈ Z, we know that f∗j = LrKj ,h(g ∗ j,h) for some g ∗ h,j =∑
i≥1 βh,iψ j h,i ∈ L2(ρ T+1 Xj ). Then we have
f∗j = L r Kj ,T+1( ∑ i≥1 βh,iψ j h,i) = ∑ i≥1 (ζjh,i) rβh,iψ j h,i.
The case r = 12 means each f ∗ j lies in the RKHSHKj . Then we have
f̃η,h,j − f∗j = (LKj ,h + ητjI)−1LKj ,hf∗j − f∗j = ∑ i≥1 ητj(ζ j h,i) rβjh,iψ j h,i ζji + ητj .
Then we have ‖f̃η,T+1,j − f∗j ‖22 = ( ∑ i≥1 ητj(ζ j h,i) rβjh,i ζjh,i + ητj )2 = (ητj) 2r ∑ i≥1 ( ητj ζjh,i + ητj )2−2r( ζjh,i ζjh,i + ητj )2r(βjh,i) 2
≤ (ητj)2r ∑ i≥1 (βjh,i) 2 = (ητj) 2r‖L−rKj ,T+1f ∗ j ‖22.
Similarly, we also have ητj‖f̃η,T+1,j‖22 = (ητj)2 ∑ i≥1 ( (ζjh,i) 1+r ζjh,i + ητj )2(βjh,i) 2 = (ητj) 2r ∑ i≥1 ( ζjh,i ζjh,i + ητj )2+2r( ητj ζjh,i + ητj )2−2r(βjh,i) 2
≤ (ητj)2r‖L−rKj ,T+1f ∗ j ‖22.
Proposition 3. Under Assumptions 3, there holds
E2 = E(f∗η,T+1)− E(f∗) + η p∑ j=1 τj‖f∗η,T+1,j‖2K ≤ (p+ 1) p∑ j=1 (ητj) 2r‖L−rKj ,T+1f ∗ j ‖22.
Proof. According to Lemma 1, we have
ET+1(f∗η,T+1)− ET+1(f∗) + η p∑ j=1 τj‖f∗η,T+1,j‖2K
≤ EσT+1(f∗η,T+1)− EσT+1(f∗) + η p∑ j=1 τj‖f∗η,T+1,j‖2K + Mc σc
≤ EσT+1(f̃η,T+1)− EσT+1(f∗) + η p∑ j=1 τj‖f̃η,T+1,j‖2K + Mc σc
≤ ET+1(f̃η,T+1)− ET+1(f∗) + η p∑ j=1 τj‖f̃η,T+1,j‖2K + 2Mc σc
Then we have ET+1(f̃η,T+1)− ET+1(f∗) + η p∑ j=1 τj‖f̃η,T+1,j‖2K = ‖f̃η,T+1 − f∗‖22 + η p∑ j=1 τj‖f̃η,T+1,j‖2K
≤ ‖ p∑ j=1 [f̃η,h,j − f∗j ]‖22 + η p∑ j=1 τj‖f̃η,h,j‖2K
≤ p∑ j=1 p‖f̃η,T+1,j − f∗j ‖22 + η p∑ j=1 τj‖f̃η,T+1,j‖2K
≤ (p+ 1) p∑ j=1 (ητj) 2r‖L−rKj ,T+1f ∗ j ‖22.
Combing above results, we get the desired result.

heading: B.3 HYPOTHESIS ERROR
text: This section focuses on bounding the hypothesis error
E3 = {EσT (f̂) + λΩ(f̂)− EσT (f̂η)− η p∑ j=1 ‖f̂η,j‖2K}.
We first give the properties of f̂ and then use them to bridge f̂ and fη . Lemma 6. For all j = 1, ..., p, there holds
τj‖αηj ‖2 = 1
ηT √√√√ T∑ t=1 (`′σ(y t − p∑ j=1 K ′jtα η j )) 2 ≤ 2σ ηT 1 2
Proof. Recall the represent theorem which ensures that
f̂η = p∑ j=1 T∑ t=1 αηtjKj(x t j , ·), α η tj ∈ R.
For notation simplicity, denote αηj = (α η 1j , ..., α η Tj) ′ ∈ RT and αη = ((αη1)′, ..., (αηp)′)′ ∈ RTp. From the definition of Eq. (7), we deduce that
αη = arg min α∈RTp { 1 T T∑ t=1 `σ(y t − p∑ j=1 K ′jtαj) + η p∑ j=1 τj(α η j ) ′Kjαj},
whereKj = {K(xsj , xij)}ni,s=1 ∈ RT×T . we then have
1
T T∑ t=1 `′σ(y t − p∑ j=1 K ′jtα η j )K ′ ji = ητjKjα η j .
It is easy to deduce that
αηj = 1 ητjT (`′σ(y 1 − p∑ j=1 K ′j1α η j ), ..., ` ′ σ(y T − p∑ j=1 K ′jTα η j )) ′.
Then we obtain that, for any j = 1, ..., p,
τj‖αηj ‖2 = 1
ηT √√√√ T∑ t=1 (`′σ(y t − p∑ j=1 K ′jtα η j )) 2 ≤ 2σ ηT 1 2
This completes the proof.
Proposition 4. Under Assumptions 3, there holds
E3 ≤ λΩ(f̂η) = λ p∑ j=1 τj‖αηj ‖2 ≤ 2λpσ ηT 1 2 .

heading: Proof.
text: E3 = {EσT (f̂) + λΩ(f̂)− EσT (f̂η)− η p∑ j=1 ‖f̂η,j‖2K}
= EσT (f̂) + λΩ(f̂)− EσT (f̂η)− λΩ(f̂η)− η p∑ j=1 ‖f̂η,j‖2K + λΩ(f̂η)
≤ λΩ(f̂η) Combining the above inequality with Lemma 6, we get that
E3 ≤ λΩ(f̂η) = λ p∑ j=1 τj‖αηj ‖2 ≤ 2λpσ ηT 1 2 .
Proof of Theorem 1: Combining Propositions 1-4, for Tη ≤ 1, we have with confidence 1− δ
‖f−f∗‖2 L2(ρ T+1 X ) ≤ C̃ log(1/δ)(T− 12 +T−1σ+η2r+λη−1T− 12σ+σ−cλ−2T−2 +σ−cλ−1T−1),
where C̃ is a positive constant independently of T, λ, η, δ, σ. By taking σ = Tm, η = T β and λ = T γ , we then have with confidence 1− δ
‖f − f∗‖2 L2(ρ T+1 X ) ≤ C̃ log(1/δ)TΨ(m,β,γ,c,r),
where
Ψ(m,β, γ, c, r) = max{−1 2 ,m− 1, 2rβ, γ − 1 2 − β +m,−cm− 2γ − 2,−cm− γ − 1}.
By taking σ = 12c , β = − 1 4r and γ = − 1 4r − 1 2c Direct computation shows that
Ψ(m, c, r) =
{ max{− 12 , 1 2c − 1}, if m ≤ 1− 1 4r
max{− 12 ,m− 1,−cm+ 1 2r + 2m− 2}, if m > 1− 1 4r .
This completes the proof.

heading: C PROOF OF THEOREM 2
text: Denote by s = {st}Tt=1 a probability set with ∑T t=1 st = 1. We define following functions associated with the probability set s:
f∗T+1 = p∑ j=1 f∗T+1,j = arg min f∈HK E`σ(f(xT+1)− yT+1),
f∗,sT = p∑ j=1 f∗,sT,j = arg min f∈HK T∑ t=1 st`σ(f(x t)− yt),
f̂ sη = p∑ j=1 f̂ sη,j = arg min f= ∑p j=1 fj ,fj∈HKj { T∑ t=1 st`σ(y t − p∑ j=1 fj(x t j)) + η p∑ j=1 τj‖fj‖2Kj},
and
f̂ s = p∑ j=1 f̂ sj = arg min f= ∑p j=1 fj ,fj∈HKj { T∑ t=1 st`σ(y t − p∑ j=1 fj(x t j)) + λΩ(f)}.
Correspondingly, we denote by Eσ,sT (f) = ∑T t=1 st`σ(f(x t)− yt).
Proposition 5. Let Assumptions 1-2 be true. For any f ∈ H̄Z , there holds
‖f − f∗‖2 L2(ρ T+1 X ) ≤ E1 + E2 + E3 + EσT (f∗T ) + η p∑ j=1 τj‖f∗T,j‖2Kj + 2Mc σc ,
where E1 = EσT+1(f̂ s)− E σ,s T (f̂ s)
E2 = Eσ,sT (f̂ s) + λΩ(f̂ s)− Eσ,sT (f̂ s η)− η p∑ j=1 τj‖f̂ sη,j‖2Kj
and
E3 = ET+1(f∗η,T+1)− ET+1(f∗) + η p∑ j=1 τj‖f∗η,T+1,j‖2Kj .
Proof. According to Lemma 2, for any f ∈ H̄Z , we can make following error decomposition
‖f − f∗‖2 L2(ρ T+1 X )
= ET+1(f)− ET+1(f∗)
≤ EσT+1(f)− EσT+1(f∗) + Mc σc ≤ {EσT+1(f)− Eσ,sT (f)}︸ ︷︷ ︸ E1 +{Eσ,sT (f)− E σ,s T (f ∗,s T )}+ E σ,s T (f ∗,s T ) + E σ T+1(f ∗ η,T+1)− EσT+1(f∗)
+ η p∑ j=1 τj‖f∗η,T+1,j‖2Kj}+ Mc σc
≤ E1 + {Eσ,sT (f) + λΩ(f)− E σ,s T (f ∗,s T )− η p∑ j=1 τj‖f∗,sT,j‖ 2 Kj}+ E σ,s T (f ∗,s T ) + η p∑ j=1 τj‖f∗,sT,j‖ 2 Kj
+ EσT+1(f∗η,T+1)− EσT+1(f∗) + η p∑ j=1 τj‖f∗η,T+1,j‖2Kj}+ Mc σc
≤ E1 + {Eσ,sT (f) + λΩ(f)− E σ,s T (f̂ s η)− η p∑ j=1
τj‖f̂sη,j‖2Kj}︸ ︷︷ ︸ E2 +Eσ,sT (f ∗,s T ) + η
p∑ j=1 τj‖f∗,sT,j‖ 2 Kj
+ EσT+1(f∗η,T+1)− EσT+1(f∗) + η p∑ j=1 τj‖f∗η,T+1,j‖2Kj}+ Mc σc
≤ E1 + E2 + ET+1(f∗η,T+1)− ET+1(f∗) + η p∑ j=1
τj‖f∗η,T+1,j‖2Kj︸ ︷︷ ︸ E3 +Eσ,sT (f ∗,s T ) + η
p∑ j=1 τj‖f∗,sT,j‖ 2 Kj + 2Mc σc
≤ E1 + E2 + E3 + Eσ,sT (f ∗,s T ) + η p∑ j=1 τj‖f∗,sT,j‖ 2 Kj + 2Mc σc
≤ E1 + E2 + E3 + min f∈HZ {Eσ,sT (f) + λΩ(f)}+ η p∑ j=1 τj‖f∗,sT,j‖ 2 Kj + 2Mc σc .
Following the generalization analysis in (Kuznetsov & Mohri, 2020a), we need to discrepancy measure to measure the discrepancy of the non-stationarity of the stochastic process {Zt}∞t=−∞ with respect to both the loss function `σ and the hypothesis set H̄Z .
Definition 4. The discrepancy describes the discrepancy between target distribution and the distribution of the sample. For any f ∈ H̄Z , the discrepancy measure with respect to Huber loss is defined as
disc(s) := sup f∈HZ
{ E`σ(f(xT+1)− yT+1)− T∑ t=1 Est`σ(f(xt)− yt) } .
Lemma 7. For any δ > 0, with probability at least 1 − δ, the following inequality holds for all f ∈ H̄Z and all α > 0 E`σ(f(xT+1)−yT+1) ≤ T∑ t=1 st`σ(f(x t)−yt)+disc(s)+‖s‖2+6M`σ √ 4π log TRT (G)+M`σ‖s‖2 √ 8 log 1 δ ,
where G := {`σ(f, z) : f ∈ H̄Z} and M`σ is a positive constant such that `σ(f(xt)− yt) ≤ M`σ for any f ∈ H̄Z and t ∈ Z.
Proof of Theorem 2 The proofs of bounding errors E1, E2 and E3 proceeds similarly to the proof of Theorem 1 and are omitted for brevity.
According to Lemma 3, we have
E1 ≤ disc(s) + ‖s‖2 + 6M`σ √ 4π log TR+M`σ‖s‖2 √ 8 log 1
δ ,
where R = ‖s‖2[32σ[1 + 4 √ 2 log3/2(eT 2)]M2p 1 2κ].
Moreover, from Proposition 4, we can get
E3 ≤ (p+ 1) p∑ j=1 (ητj) 2r‖L−rKj ,T+1f ∗ j ‖22.
Similarly, to bound the hypothesis error
E2 = {Eσ,sT (f̂ s) + λΩ(f̂ s)− Eσ,sT (f̂ s η)− η p∑ j=1 ‖f̂sη,j‖2K},
we first give the properties of f̂ and then use them to bridge f̂ and fη .
Lemma 8. For all j = 1, ..., p, there holds
E2 ≤ 2λpσ‖s‖2
η .
Proof. Recall the represent theorem which ensures that
f̂ sη = p∑ j=1 T∑ t=1 αη,stj Kj(x t j , ·), α η,s tj ∈ R.
For notation simplicity, denote αη,sj = (α η,s 1j , ..., α η,s Tj ) ′ ∈ RT and αη,s = ((αη,s1 )′, ..., (αη,sp )′)′ ∈ RTp. We then deduce that
αη,s = arg min α∈RTp { T∑ t=1 st`σ(y t − p∑ j=1 K ′jtαj) + η p∑ j=1 τj(α η j ) ′Kjαj},
whereKj = {K(xsj , xij)}ni,s=1 ∈ RT×T . We get
T∑ t=1 st` ′ σ(y t − p∑ j=1 K ′jtα η,s j )K ′ ji = ητjKjα η,s j .
It is easy to deduce that
αη,sj = 1
ητj (s1`
′ σ(y 1 − p∑ j=1 K ′j1α η j ), ..., sT ` ′ σ(y T − p∑ j=1 K ′jTα η j )) ′.
Then we obtain that, for any j = 1, ..., p,
τj‖αηj ‖2 = 1
η √√√√ T∑ t=1 s2t (` ′ σ(y t − p∑ j=1 K ′jtα η j )) 2 ≤ 2σ‖s‖2 η
This completes the proof.
Under Assumptions 1-2, there holds
E2 = {EσT (f) + λΩ(f̂ s)− EσT (f̂ sη)− η p∑ j=1 ‖f̂ sη,j‖2K}
= EσT (f̂ s) + λΩ(f̂ s)− EσT (f̂ sη)− λΩ(f̂ sη)− η p∑ j=1 ‖f̂ sη,j‖2K + λΩ(f̂ sη)
≤ λΩ(f̂ sη) = λ p∑ j=1 τj‖αη,sj ‖2
≤ 2λpσ‖s‖2 η .
By combining the above results, we have with confidence 1− δ
‖f̂ s − f∗‖2 L2(ρ T+1 X )
≤ disc(s) + min f∈HZ {Eσ,sT (f) + λΩ(f)}
+C̃ log(1/δ)(‖s‖2 + ‖s‖2σ + η2r + λη−1‖s‖2σ + σ−cλ−2T−2 + σ−cλ−1T−1 + η),
where C̃ is a positive constant independently of T, λ, η, δ, σ and p. By taking σ = T 1 2c , λ = T−1 and η = T− 1 2 , we completes the proof.

heading: D PROOF OF THEOREM 3
text: According to the definition of disc(T + 1), we have
disc(s) = sup f∈HZ
{ E`σ(f(xT+1)− yT+1)− T∑ t=1 stE`σ(f(xt)− yt) }
≤ sup f∈HZ [ T∑ t=1 q∗tE`σ(f(xt)− yt)− T∑ t=1 stE`σ(f(xt)− yt)] + disc(q∗).
Furthermore, it is easy to deduce that
sup f∈HZ T∑ t=1 (q∗t − st)E`σ(f(xt)− yt)− sup f∈HZ T∑ t=1 (q∗t − st)`σ(f(xt)− yt)
≤ sup f∈HZ T∑ t=1 (q∗t − st)[E`σ(f(xt)− yt)− `σ(f(xt)− yt)].
Proof of Theorem 3 According to the proof of Theorem 1 in Kuznetsov & Mohri (2020a), we have
sup f∈HZ T∑ t=1 (q∗t−st)[E`σ(f(xt)−yt)−`σ(f(xt)−yt)] ≤ 1√ T +6M`σ √ 4π log TRT (G)+M`σ‖q∗−s‖2 √ 8 log 1 δ
By combining above results, we have with confidence 1− δ ‖f̂ s − f∗‖2
L2(ρ T+1 X )
≤ disc(q∗) + sup f∈HZ T∑ t=1 (q∗t − st)`σ(f(xt)− yt) + min f∈HZ {Eσ,sT (f) + λΩ(f)}
+C̃ log(1/δ)(‖q∗ − s‖2 + ‖s‖2 + ‖s‖2σ + η2r + λη−1‖s‖2σ + σ−cλ−2T−2 + σ−cλ−1T−1 + η).
By taking σ = T 1 2c , λ = T−1 and η = T− 1 2 , we complete the proof.
Algorithm 1: Optimization procedure for adaptive SpHAM Input: Data {(xt, yt)}Tt=1, Max-Iter Z ∈ Z, Mercer kernel Kj , j = 1, ..., p with bandwidth d, Weights τl, l = 1, ..., p, q∗. Initialization: Lipschitz constant L, s0. Step A: Computing weights ŝ: for z = 1, ..., Z do
1. Compute Az−1 via DC-programming; 2. Update sz via (12).

heading: Output: ŝ = sZ . Step B: Computing f̂ s:
text: for z = 1, ..., Z do
1): Compute αz = pL(βz) via (14);
2): mz+1 = 1+ √ 1+4m2z 2 ; 3): βz+1 = αz + mz−1mz+1 (α z − αz−1).
Output: αŝ = αZ ; Prediction function: f̂ s = ∑P j=1 ∑T t=1 α ŝ] tjKj(x t j , ·); Variable selection: {j : ‖αsj‖2 ≥ v, j = 1, ..., p}.

heading: E EXPERIMENT OPTIMIZATION
text: The optimization problem (10) reduces to problem 7 when taking s∗ = 1T IT . Recall the optimization problem in Step A:
ŝ = arg min s { sup f∈HZ T∑ t=1 (q∗t − st)`σ(f(xt)− yt) + λ1 2 ‖q∗ − s‖22 + λ2 2 ‖s‖22T 1 2c }.
The above optimization problem can be equivalently rewritten as a common type of bilevel optimization problem (Colson et al., 2007), i.e,
Outer problem: mins ∑T t=1(q ∗ t − st)`σ(fs(xt)− yt) + λ1‖q∗ − s‖22 + λ2T 1/2c‖s‖22,
Inner problem: fs = arg maxf∈Hz ∑T t=1(q ∗ t − st)`σ(f(xt)− yt),
where the outer (min) problem parameterized by s, is nested within the inner (max) problem. The outer problem can be solved by standard gradient descent, where in each step, we need to optimize the inner (max) problem with last updated s.
We denote by k the iteration time. Then for k + 1-th iteration, we have the following gradient update rule
sk+1 = sk − γ(Ak − λ1(q∗ − sk) + λ2T 1 2c sk), (12)
where γ is learning rate,
Ak = − ( `σ( T∑ t=1 p∑ j=1 αs k tjKj(x t j , x 1 j )− y1), ..., `σ( T∑ t=1 p∑ j=1 αs k tjKj(x t j , x T j )− yT ) )′ ∈ RT
and αs k = (αs k
11, ..., α sk T1, ..., α sk 1p, ..., α sk Tp) ′ ∈ RTp
is obtained by the following weighted optimization problem
αs k
= arg max α T∑ t=1 (q∗t − skt )`σ( T∑ m=1 p∑ j=1 αmjKj(x m j , x t j)− yt).
Note that the inner problem is subjected to f ∈ HK , i.e., ‖α‖2 = ∑p j=1 τj‖αj‖2 ≤
2M2
λT minj=1,...,p τj .
A widely-used method for solving this constrained inner problem is DC programming (Tao & An, 1998). After obtaining the solution ŝ, we turn to solve the following weighted optimization problem in Step B:
αŝ = arg min αj∈RT ,j=1,...,p { T∑ t=1 ŝt`σ(y t − p∑ j=1 (Ktj) ′αj) + λ p∑ j=1 τj ||αj ||2}. (13)
We can see that this problem contains non-smooth function λ ∑p j=1 τj ||αj ||2 which makes the standard gradient descent method inapplicable. To conquer this challenge, we leverage fast iterative shrinkage-thresholding algorithm (FISTA)(Beck & Teboulle, 2009). Our optimization problem becomes
αk = pL(α k−1) := arg min
α
{ λ p∑ j τj ||αj ||2 + L 2 ||α− (αk−1 − 1 L T∑ t=1 ŝt∇`σ(yt − (Kt)′αk−1)||2 }
= (( 1− λ ||vj ||2 ) + vj ) 1≤j≤p , (14)
where
vj = α k−1 j −
1
L T∑ t=1 ŝt∇αj `σ(yt − (Kt)′αk−1))
and
∇αj `σ(yt − (Kt)′αk−1)) =  2((Kt)′α− yt)(Ktj)′ |(Kt)′α− yt| < σ 2σ (Kt)′α− yt ≥ σ −2σ (Kt)′α− yt ≤ −σ
and L = max 1T ∑T t=1 ‖(Kt)′Kt‖ is the lipschitz constant 1 T ∑T t=1∇`σ(yt − (Kt)′αk−1). Denote Proj(s) as a projection of s and v > 0 as the threshold value. To focus on the weights of samples which are really useful for forecasting, we consider Proj(s) as a box projection such that st > q∗t ,∀t = 1, ..., p). Finally, the optimization procedure for adaptive SpHAM can be summarized in Algorithm 1. Note that if we only run Step B with ŝt = 1T ,∀t = 1, ..., T , we can further obtain the optimization procedure for SpHAM. Remark 3. The computational complexity of Algorithm 1 depends on the optimization strategy for DC programming, the training size T , dimension p and the iteration times Z. We denote by the computational complexity of DC programming O(DC(T, p)). Then, the computational complexity of Step A is O(ZDC(T, p) + ZT ) and the computational complexity of Step B is O(ZT 3p3). Thus, the total computational complexity of Algorithm 1 is O(ZO(DC) + ZT 3p3). For large scale data, we can further speed up Algorithm 1 by random Fourier features technique [Rahimi and Recht 2007], which is leaved for future work. This has been carefully discussed in the revised manuscript.

heading: F ADDITIONAL EXPERIMENT

heading: F.1 EVALUATION ON HYPER-PARAMETERS
text: We investigate the impact of the selection of σ on our method. For fixed sample size T = 400, the empirical results over different σ = Tα are summarized in Table 4. This table shows that, for heavy-tailed t noise, our method tends to perform better with smaller α. This result also validates Theorem 1 and Figure 1.

heading: F.2 EVALUATION ON BENCHMARK DATA
text: We test our algorithm on nonlinear dataset from CauseMe. The hyper-parameter selection is the same as the one in Section 4. The results in Table 3 verify the effectiveness of our method.

heading: F.3 EXPERIMENTS ON AIR QUALITY DATASET
text: We use the Air Quality dataset obtained from UCI Machine Learning Repository (https:// archive.ics.uci.edu/ml/datasets/Air+quality) to test our model’s ability to detect the Granger causality. This dataset includes 9358 hourly air quality data in an Italian city, collected from March 2004 to February 2005. The details of the dataset can be obtained on the UCI website. The Granger causal network we detect is shown in Figure 3. Given the network, we observe that the temperature (T) influences ozone (O3) and nitrogen dioxide (NO2) which is validated in Kalisa et al. (2018). Mwaniki et al. (2014) confirm the relationship between relative humidity (RH) and nitrogen dioxide (NO2). Yan et al. (2018) verify the relationship between humidity and nitrogen dioxide (NO2).

heading: F.4 EXPERIMENTS ON CORONAL MASS EJECTIONS DATASET
text: Coronal Mass Ejections (CMEs) are the most violent eruptions in the Solar System. Despite machine learning approaches have been applied to these tasks recently Wang et al. (2019); Liu et al. (2018), there is no any work for interpretable prediction with Granger causal network. CMEs data are provided in The Richardson and Cane List (http://www.srl.caltech.edu/ACE/ ASC/DATA/level3/icmetable2.htm). From this link, we collect 152 ICMEs observations from 1996 to 2016. The features of CMEs are provided in SOHO LASCO CME Catalog (https: //cdaw.gsfc.nasa.gov/CME_list/). In-situ solar wind parameters can be downloaded from OMNIWeb Plus (https://omniweb.gsfc.nasa.gov/). A total of 9 features are chosen as input, including: (1) Central PA (CPA), (2) Angular Width, (3) three approximated speeds ( Linear Speed, 2nd-order Speed at final height,and 2nd-order Speed at 20 Rs), (4) Mass, (5) Kinetic Energy, (6) MPA and (7)CMEs arrival time. Figure 4 shows that Granger causal network when the output is CMEs Arrival time. Some interesting findings are concluded from this Granger causal network. For instance, Speed and Mass, as the significant variables causing CMEs arrival time, have been also screened out in Liu et al. (2018). Morover, the CMEs Angular Width does not cause the CMEs arrival time forecasting, while Liu et al. (2018) state that they have a significant correlation. This indicates that the CMEs arrival time is affected by the Angular Width at current time, but not by the historical Angular Width.

Review summaries:
This ms tackles a timer series problem where each observation is the sum of p hidden functions and a noise (termed as innovation in ms). Traditionally the noise is Gaussian, but this ms makes the noise into heavy tailed distribution such as t distribution. The L2 loss function of the original problem is then replaced by so called Huber loss, which introduces a cutoff to the deviations and is less affected by outliers than L2 loss. Based on this new loss function, this ms sets up the sparse Huber additive model and derive the bounds of the learned hidden functions, under stationary processes. In case of non-stationary processes, the ms introduces different weights to each time points, which can be used to calculate a discrepancy that reflects the non-stationarity. With different weights to different time points, the ms come up with new error bounds. Inference algorithms are provided in the end.

The authors proposed a robust Huber additive model for non-stationary time series prediction. They combine the idea of robust Huber regression (against non-Gaussian innovations) with a linear additive model formulated with the aid of the representer theorem in RKHS; sparsity is imposed on the weights, and the overall weights optimization problem can be solved efficiently through the classic FISTA algorithm due to the special problem structure. The major contributions lie in a robust (due to the Huber loss function) and interpretable (due to the use of linear additive model) time series prediction model as well as some theoretical supports on the performance. 

The authors propose, theoretically analyze and empirically evaluate a sparse additive huber additive model (SpHAM), which is an additive model where each component function in the model is a kernel function that is fitted by minimizing a huber huber loss function with a functional sparsity inducing penalty term.  Statistical properties when estimating this model on stationary data with non gaussian noise and also on non-stationary data are developed. Finally, the model and fitting routines are demonstrated empirically on 2 synthetic and 1 real world time series.