A sparsity-controlled vector autoregressive model

A Sparsity-Controlled Vector Autoregressive Model Emilio Carrizosaa, Alba V. Olivares-Nadal 1,a and Pepa Ram´ırez-Cobob aInstituto de Matem´aticas de la Universidad de Sevilla (Spain), bDepartment of Statistics and Operational Research, Universidad de C´adiz (Spain), Abstract Vector autoregressive (VAR) models constitute a powerful and well studied tool to analyze multivariate time series. Since sparseness, crucial to identify and visualize joint dependencies and relevant causalities, is not expected to happen in the standard VAR model, several sparse variants have been introduced in the literature. However, in some cases it might be of interest to control some dimensions of the sparsity, as e.g. the number of causal features allowed in the prediction. To authors extent none of the existent methods endows the user with full control over the different aspects of the sparsity of the solution. In this paper we propose a sparsity-controlled VAR model which allows to control different dimensions of the sparsity, enabling a proper visualization of potential causalities and dependencies. The model coefficients are found as the solution to a mathematical optimization problem, solvable by standard numerical optimization routines. The tests performed on both simulated and real-life multivariate time series show that our approach may outperform both the standard and Group Lasso in terms of prediction errors specially when highly sparse graphs are sought, while avoiding the VAR’s overfitting for more dense graphs. Causality; Mixed Integer Non Linear Programming; multivariate time series; sparse models; Vector autoregressive process. 1Corresponding author: Alba V. Olivares-Nadal, Departamento de Estad´ıstica e Investigaci´on Operativa, Facultad de Matem´aticas, Universidad de Sevilla, Av. Reina Mercedes, s/n, 41012 Sevilla, Spain. Phone number: +34 637 872 829 Email address: [email protected] (Alba V. Olivares-Nadal) 1 1 Introduction A plethora of real world data, such as e.g. air pollution measures or brain functional connections, involve multivariate time series, i.e., different and inter-related features, evolving in time, are simultaneously measured and are to be forecasted. Since the components of such multivariate time series are not independent, inaccurate predictions are expected if the series are analyzed separately by repeatedly using one-dimensional time series forecasting tools. In order to properly address such dependencies, vector autoregressive models (VAR) are frequently applied. However, although capturing features dependencies, there is no reason to expect the so-obtained VAR to be sparse. In other words, the output may be too complex when, on top of obtaining sharp forecasts, visualization of relevant causalities is sought, (Eichler, 2012; Shojaie and Michailidis, 2010; Lozano et al., 2009; Arnold et al., 2007; Vald´es-Sosa et al., 2005). Visualization is a trending topic nowadays, and it is being widely studied to unravel potential causalities in biological systems. In particular, graphical models have been developed to deal with genetic networks (Abegaz and Wit, 2013; Shojaie and Michailidis, 2010; Lozano et al., 2009; Hu and Hu, 2009; Dobra et al., 2004), which includes E. coli and Arabidopsis thaliana regulatory networks, or human cancer cell data. Moreover, Gorrostieta et al. (2013); Vald´es-Sosa et al. (2005) develop sparse autoregressive models to enhance visualization of brain functional connectivity. Also, understanding relevant causalities has played an important role in evaluating the effect of air pollution and exposure over human health (Dominici et al., 2000). Due to its wide range of applications, several attempts have been proposed in the literature to obtain VAR models in which sparsity, as a potential for easy visualization of complex causal relations, is pursued. This is the case, for instance, of Stochastic Search Variable Selection (George, 2000; George and McCulloch, 1997), Bayesian approaches (Doan et al., 1984) and the Lasso (Tibshirani, 1996). Minimizing the forecasting errors plus an `1-penalty regularization term has not led to a unique Lasso method, but it has evolved into a full class of Lasso approaches, such as the so-called Adaptative Lasso, (Zou, 2006), the Group Lasso, (Song and Bickel, 2011; Haufe et al., 2010; Zhao et al., 2009; Yuan and Lin, 2006), the Maximum Likelihood Estimated Lasso (Davis et al., 2012; Hsu et al., 2008), or Lasso Granger methods, (Arnold et al., 2007). Such VAR models attempt to gain overall sparsity, without an explicit analysis of its different aspects. Nevertheless, it might be of interest in certain real life situations to control not only the overall number of depicted dependencies but also other levels of sparsity. For example, the number of causal features might be wanted to be limited when acceding to the historical records of the features incurs into a cost. Consider a patient who is under surveillance and several tests need to be undertaken to control 2 her health periodically. Those tests may not only be costly but also invasive for the patient end therefore it may be desirable to reduce their application without affecting much the quality of the diagnose. As an example, Griffin et al. (2005) suggest heart rate to predict neonatal sepsis, instead of obtaining blood from the infants for laboratory tests. This is also useful for the companies shares prices, whose time series are to be paid but, once done, all the historical records are available to use. In all these cases the number of causal features and the number of dependencies from each causal feature are valued differently: as we might seek to limit the first one, we might not be that strict with the second. This subject was previously noted by Lozano et al. (2009), who stated that “as a method of Granger graphical modeling, the relevant variable selection question is not wether an individual lagged variable is to be included in regression, but wether the lagged variables for a given time series as a group (i.e. the feature), are to be included”. To address this issue they proposed to use Group Lasso, in which all the lagged variables of a feature were assigned to the same group. Although this approach reinforces to choose all the past values of a feature once one of them has already been selected for the prediction, it still does not grant control over the exact level of sparsity of the outcome, like none of the other Lasso methods. Moreover, the Lasso approaches perform a shrinkage over the coefficients which, as will be seen in Section 4, might not be advisable when highly sparse graphs are sought. In contrast we propose a novel sparse approach, formulated as a mathematical optimization problem, which endows the user with the power to control different aspects of the sparsity of the solution. Sparsity is meant here in some of its many different dimensions, as the total number of nonzero coefficients, the total number of features used for the forecast, or the number of past observations of each feature used by the model to make predictions. The performance of this sparsity-controlled VAR method (SC-VAR henceforth) will be compared with three benchmark approaches, namely, the VAR and both the standard and Group Lasso, on simulated and real-life multivariate time series. The results show that the proposed approach outperforms the benchmark Lasso methods in terms of prediction errors when highly sparse graphs are sought, while avoiding the VAR’s overfitting for more dense graphs. The paper is structured as follows. Next section introduces mathematically the three benchmark methods, the VAR, the classic Lasso and the Group Lasso, and motivates our approach. In Section 3 the SC-VAR model is introduced and expressed as a mixed integer non-linear optimization program, solvable by standard optimization software. Also, a discussion about the choice of the parameters of the model is included. Competing approaches are compared against the proposed method in Section 4 in both simulated and real datasets. Finally, conclusions and 3 future extensions are collected in Section 5. 2 Preliminaries Let {Xt}t≥0be an N-dimensional vector autoregressive process of order p, VAR(p), i.e., each series i,i= 1, ..., N, can be expressed as Xi,t =ci+ N X j=1 p X k=1 αi jkXj,t−k+ei tt≥0 where {ei t}t≥0denotes the series of contemporaneous shocks that affect feature i, and ciand αi jk are real numbers. The usual estimation procedures for the coefficients ciand αi jk are Maximum Likelihood (which implies making distributional assumptions on the errors ei t) or, without imposing any statistical assumption, the Ordinary Least Squares method: min c,A N X i=1 T X t=p Xi,t+1 −ci− N X j=1 p X k=1 αi jkXj,t+1−k!2 (1) where c= (ci)i∈R1×Nand A= (A1|...|AN)∈RN×Np stand for all unknown coefficients of the process to be estimated. Here Ai= (αi jk)j,k represents the N×p matrix of coefficients used to model series i. There is no reason to expect sparsity in the estimates obtained by maximum likelihood estimation or by solving the nonlinear program (1), and therefore, it may be difficult to visualize causalities while leading to overfitting (Kalli and Griffin, 2014; Li, 2012; Kojima et al., 2009). Among the different procedures proposed in the literature with the aim of obtaining more sparse solutions, a prominent role is given to the Lasso-VAR (Lasso thereafter), in which an `1regularization term is added to the objective function (1), and thus estimates are obtained by solving the following nonlinear nonsmooth optimization program: min c,A N X i=1 T X t=p Xi,t+1 −ci− N X j=1 p X k=1 αi jkXj,t+1−k!2 + N X i=1 λi N X j=1 p X k=1 |αi jk|!.(2) As previously mentioned, when features are presented as time series, the use of Group Lasso is encouraged in the literature if more than one-lagged values are considered; see, for instance, Lozano et al. (2009). Such a method groups lagged variables of the same feature, giving more importance to the number of causal features rather than to the overall number of dependencies depicted. This is done by adding an `2-penalty separately for each group (Yuan and Lin, 2006): 4 min c,A N X i=1 T X t=p Xi,t+1 −ci− N X j=1 p X k=1 αi jkXj,t+1−k!2 + N X i=1 λi R X r=1 kαi j(r) αi j(r) αi j(r)k2!(3) where αi j(r) αi j(r) αi j(r)=αi j1, αi j2, ..., αi jp. Infinitely many solutions may be obtained for the Lasso approaches when varying the parameter λi∈R+, where more sparse solutions are attained when increasing the value of this penalty (in fact, note that all αi jk = 0 when λi→+∞). However, although one may think that these popular Lasso approaches already cope with sparsity-controlled graphs, there is no way to discern the exact level of sparsity of the outcome when values for the penalties λiare given a priori, calling for a (multidimensional) parameter tuning which, at the end, show solutions with different levels of sparsity, but there is no guarantee that sparsity is fulfilled. Moreover, these methods perform a shrinkage over the coefficients, which may deteriorate their prediction power when highly sparse graphs are sought. Even although these methods often provide reasonable sparse solutions, according to some authors (Bertsimas and Copenhaver, 2014; Caramanis et al., 2012) adding an `1or `2penalty is a robust approach rather than a sparse method. In contrast, Bertsimas and Copenhaver (2014) present Mixed Integer Programming as a useful tool to attain sparsity. This recommendation is supported by the recent improvement on integer optimization solvers, which can attain considerably good solutions at a reasonable computational cost. In this paper a sparsity controlled VAR is formulated as a Mixed Integer Non-Linear Problem (MINLP), in which we manage to control different aspects of the sparsity of the outcome. When control over sparsity is sought, a natural approach is as follows: start with a standard VAR and then, in a naive way, select the largest estimated coefficients (in absolute value) and set to zero the remaining (smaller) coefficients. This naive approach is easy to implement, quick to execute and sounds reasonable. However, its performance may be poor, as we show next. Consider the top left panel of Figure 1, where a simulated VAR of order p= 3 is graphically represented. The multivariate time series is represented as a directed graph; nodes correspond to the different features, i.e., the different one-dimensional time series composing the multivariate time series; edges in the graph visualize causality: an arrow from node jto node imeans that the model uses feature jin the forecast of feature i; the thickness of the edges is proportional to the magnitude (in absolute value) of the coefficient relating the features, and therefore it measures the causality’s strength when the series are normalized. The color of the edges is related to the lag: the arrow is plotted in black if the present (t) value of a feature is related with a data one period behind (t−1), red for two periods (t−2) and green for three (t−3). Here node 1 receives arrows from nodes 2 and 3, meaning that, in 5 order to forecast feature 1, past values of features 2 and 3 are used. Note also that node 1 receives here three arrows from node 3, so the present value tof feature 1 is caused by the previous three values (t−1), (t−2) and (t−3) of series 3. Assume that only one observation of one single feature is allowed to be used to explain a feature. Then the naive approach, illustrated in the top right panel of Figure 1, underestimates strong persistence in favor of just one significant coefficient. This undesired phenomenon is a consequence of the nature of the naive approach: arrows are considered to be kept or removed one by one, and thus the overall picture may be lost. For this reason, instead of using the above-described naive approach, we suggest to express the problem of arrows selection as a mathematical optimization model, solvable by current standard numerical optimization routines. In particular, when applied to the time series of the example, the output of our procedure is visualized through the directed graph in Figure 1 (bottom right panel). It can be observed that the MSE obtained by the SC-VAR solution is considerably smaller than that of the naive approach. In order to illustrate the effect of the shrinkage of Lasso methods over the forecasting power we depict the solution under the standard Lasso in the bottom left panel. This effect is clear since the arrow is thinner, which could be detrimental when highly sparse solutions are sought. This phenomenon, also observed in the results of the experiments carried out in Section 4, will be discussed in more detail. Also it is interesting to note that the chosen causal feature differs for each method: as variable 1 is considered to cause itself for the Group Lasso, feature 3 is considered for our approach instead. This is not surprising since it is known that the VAR might not be identifiable (L¨utkepohl, 2005). 3 The SC-VAR. In this section we will first discuss the sparsity parameters in the SC-VAR model (Section 3.1), which will be later written as a mathematical optimization problem (Section 3.2). In Section 3.3 we will briefly discuss the choice of parameters of our model. 3.1 Different aspects of the sparsity To the best of our knowledge the existent sparse VAR models attempt to gain overall sparsity, without an explicit analysis of the different aspects of sparsity. The Lasso approaches do not grant the user with the ability to manage the sparsity of the solution in some desirable dimensions either. Indeed, the number of features allowed to be used by the model to explain feature i(Vi S) or the overall number of nonzero coefficients (VA), are different measures usually masked under the generic 6 term of sparsity. Other aspects of the sparsity, that will be also under the user’s control in our proposed sparse VAR, are the number of non-zeroes used to explain feature i(Vi T) or the number of observations per each causal feature of variable i(Vi Sa). Moreover, when series are normalized then the strength of a potential causality might be related with the magnitude of its associated coefficient. Hence, in order to avoid spurious dependencies, clear cut-offs i jwill be introduced, so that only coefficients with an absolute value greater than or equal to i jare allowed when relating feature iwith feature j. All these parameters, included in our SC-VAR, allow the user to obtain an output with the desired level and structure of sparsity. 3.2 Mathematical Programming formulation The objective of the SC-VAR approach if twofold: control the sparsity of the solution while not damaging much the forecasting capacity. Therefore, in our SC-VAR model, the VAR estimates are obtained by solving the optimization problem (1) imposing on the coefficients of Athe bounds represented by the sparsity parameters VA, V i T, V i S, V i Sa and i jabove. They can be expressed as linear constraints by adding logical (binary) variables. Indeed, define the variables δi jk to indicate whether a coefficient αi jk is zero or not, and variables γi jto indicate if feature jis meant to cause variable i(i.e., if αi jk 6= 0 for some k). Now the SC-VAR model is formulated as the optimization problem (4)-(11), whose outputs are the estimates of the sparse coefficients ciand αi jk, as well as the 7 solution for the indicator variables δi jk and γi j. min c,A,∆,Γ, N X i=1 T X t=p Xi,t+1 −ci− N X j=1 p X k=1 αi jkXj,t+1−k!2 (4) s.t δi jk ≤γi j∀k∈K, j, i ∈I(5) N X j=1 γi j≤Vi S∀i∈I(6) p X k=1 γi jδi jk ≤Vi Sa ∀j, i ∈I(7) N X i=1 N X j=1 p X k=1 δi jk ≤VA(8) N X j=1 p X k=1 δi jk ≤Vi T∀i∈I(9) Mδi jk ≥ |αi jk| ≥ i jδi jk ∀k∈K, j, i ∈I(10) δi jk, γi j∈ {0,1} ∀k∈K, j, i ∈I(11) where K={1, ..., p},I={1, ..., N}, ∆ = (δi jk)i,j,k, Γ = (γi j)i,j and Mis a large constant. Let us briefly discuss the correctness of the formulation above. The objective function (4) minimizes the sum of squared errors. Constraint (5) forces the variable γi jto take the value 1 when some δi jk takes the value 1, i.e., as soon as some αi jk is non-zero. The remaining constraints model different aspects of the sparsity of the process. Indeed, constraints (6) and (7) bound the number of features that are said to cause variable iand the number of non-zero coefficients per each of the chosen causal features of variable i, respectively, for i= 1, .., N. Constraints (8) and (9) bound respectively the total number of non-zero entries in matrix Aand the total number of non-zero coefficients for each variable i. The shrinking parameter i jis included in the model via constraint (10), which assigns zero to any coefficient that is not allowed to appear on the model, but otherwise it requires |αi jk|to belong to the interval [i j, M]. Here Mis assumed to be a large fixed number, and thus this constraint does not exclude reasonable values of the parameters αi jk. Problem (4)-(11) is a Mixed Integer Non Linear Program (Burer and Letchford, 2012; Lee and Leyffer, 2012), called MINLP henceforth, with convex quadratic objective function. See e.g. Bertsimas and Copenhaver (2014); Bertsimas and 8 Mazumder (2014) for other statistical problems recently addressed via optimization in integer numbers. All constraints are linear, except for (10). However, this can be rewritten by introducing new auxiliary variables νi+ jk ,νi− jk , allowing to reformulate Problem (4)-(11) as: min c,A,∆,Γ, N X i=1 T X t=p Xi,t+1 −ci− N X j=1 p X k=1 αi jkXj,t+1−k!2 s.t (5)-(9), (11) αi jk ≥i jνi+ jk −(1 −νi+ jk )M∀k∈K, j, i ∈I αi jk ≤ −i jνi− jk + (1 −νi− jk )M∀k∈K, j, i ∈I δi jk ≤νi+ jk +νi− jk ≤1∀k∈K, j, i ∈I νi+ jk , νi− jk ∈ {0,1} ∀k∈K, j, i ∈I (P) Note that the new constraints require that either −M≤αi jk ≤ −i jor i j≤ αi jk ≤Mif coefficient αi jk is chosen to appear in the model (i.e., if δi jk = 1). Now, Problem (P) is a MINLP with quadratic convex objective function and linear constraints. Hence, it can be solved using standard solvers, such as CPLEX or Gurobi, and it can be easily written using a simple algebraic language such as AMPL (Fourer et al., 2002). The lines of the AMPL code are included in Appendix A of the Supplementary Material. Note that (8) is the only constraint of Problem (P) linking the count of the nonzeroes of all features. Hence, if constraint (8) were redundant, the separability of the objective function (4) would allow to solve Problem (P) by solving separately the problem for each feature i: min ci,Ai T X t=p Xi,t+1 −ci− N X j=1 p X k=1 αi jkXj,t+1−k!2 . Since such problems are of much smaller dimension, we see that removing constraint (8) would allow one to cope with databases with a larger number of time series. 3.3 Choice of paramenters Compared with other methods seeking sparsity, our approach has many more parameters, which may be, at first glance, discouraging. Some comments are in 9 bound on the total number for non-zeroes per feature (VT, left panel), and the upper-bound on the number of causal features (VS, right panel). Note that for this second case the SC-VAR is compared against the more suitable Group Lasso instead. In the left panel of Figure 6 we observe that there exists a clear difference in the impact of the parameters VSand VSa over the MSE: increasing the parameter VSa seems more efficient to reduce the prediction error than increasing VS. Therefore it is advisable to treat both parameters separately, since they clearly represent different aspects of the sparsity. This is done in a natural way with our approach, but not with the Lasso. Furthermore, note that our approach seems to stabilize for VT≥6, whatever the values of VSand VSa are. The constant MSE seems to denote that no further coefficients are being added to the model. Also, some of these constant lines explain that the constraint involving VTis redundant, as the VS·VSa ≤VT. The behaviour observed from both plots of Figure 6 is analogous to that of the previous results: when highly sparse graphs are sought, the SC-VAR seems to be more appropriate as it attains a better MSE. We point out that when only one non-zero is allowed (VT= 1), the SC-VAR yields a 16% worsening over the VAR, while the Lasso reports a 49% deterioration. Moreover, when the number of causal features wants to be limited, the differences between the prediction errors of the SC-VAR and Group Lasso approaches are roughly a 200% and 50% for one and two causal features, respectively. Note also that the MSE deterioration of the SC-VAR over the VAR is less than a 2% for VS≥2. Figure 7 depicts, through directed graphs, the solutions of the VAR, and a couple of solutions for the Group Lasso and the SC-VAR for the air pollution data set. In such graphs color blue is associated with a 4-lag dependency. From Figure 7 some comments arise. First, note that although the chosen relevant causal features are the same for the Group Lasso and SC-VAR when VS= 1, the Group Lasso implies a 191% worsening over the VAR’s prediction error, while the SC-VAR deteriorates it in a 8%. Second, the unraveled potential causal features obtained for each method are different when VS= 2. The choice of the SC-VAR provides a much more sparse graph than the VAR while obtaining 26% less deterioration than its sparse counterpart. 5 Concluding remarks and extensions In this paper a sparse vector autoregressive model, the SC-VAR, that allows the user to control the sparsity of the output from various perspectives has been introduced. The model’s sparsity is expressed in terms of different parameters, such as the number of total non-zero entries per series, the number of features involved 16 or the number of periods chosen per feature. The method is expressed as an optimization problem, solvable by standard numerical optimization software. The ability of the proposed approach to unravel potential causalities and, in many cases, to improve the fit, has been tested in simulated multivariate time series as well as in two real data bases, referred in the existent literature. It is concluded from the experiments that (i) the proposed approach is able to yield very sparse solutions either improving or without significantly increasing the VAR’s forecasting error, (ii) the SC-VAR usually yields better MSEs than both the classic and Group Lasso when highly sparse graphs are sought, leading to a much better visualization of the process dependencies, as e.g. depicted in Figure 2, and (iii) the parameters considered to measure the sparsity play different roles, thus it might not be advisable to aggregate them into one single measure, as done by other sparse methods. Although the computational times were often negligible for the simulations, where ten nodes were considered, they are expected to increase dramatically within the dimensionality of the data sets. Hence, the authors are planning to develop in the future heuristics enabling to cope with very large data sets. 6 Supplementary Material Supplementary material is available online at (to appear) Acknowledgments This research is supported by projects MTM2012-36163 (Ministerio de Econom´ıa y Competitividad, Spain), P11-FQM-7603 and FQM-329 (Junta de Andaluc´ıa), all with EU ERD Funds. We thank Richard A. Davis, Pengfei Zang and Tian Zheng for providing the concentration of air pollutants database. Conflict of Interest: None declared. References Abegaz, F. and Wit, E. (2013). Sparse time series chain graphical models for reconstructing genetic networks. Biostatistics, 14(3):586–599. Arnold, A., Liu, Y., and Abe, N. (2007). Temporal causal modeling with graphical granger methods. In Proceedings of the 13th ACM SIGKDD International Conference on Knowledge discovery and data mining, 66–75. ACM. 17 Bertsimas, D. and Copenhaver, M. S. (2014). Characterization of the equivalence of robustification and regularization in linear, median, and matrix regression. arXiv preprint arXiv:1411.6160. Bertsimas, D. and Mazumder, R. (2014). Least quantile regression via modern optimization. The Annals of Statistics, 42(6):2429–2525. Burer, S. and Letchford, A. N. (2012). Non-convex mixed-integer nonlinear programming: A survey. Surveys in Operations Research and Management Science, 17(2):97–106. Camm, J. D., Raturi, A. S., and Tsubakitani, S. (1990). Cutting big M down to size. Interfaces, 20(5):61–66. Caramanis, C., Mannor, S., and Xu, H. (2012). Robust optimization in machine learning. In Sra, S., Nowozin, S., and Wright, S. J., eds., Optimization for machine learning, 369–402. MIT Press, Michigan. Davis, R. A., Zang, P., and Zheng, T. (2012). Sparse vector autoregressive modeling. arXiv preprint arXiv:1207.0520. Doan, T., Litterman, R., and Sims, C. (1984). Forecasting and conditional projection using realistic prior distributions. Econometric reviews, 3(1):1–100. Dobra, A., Hans, C., Jones, B., Nevins, J. R., Yao, G., and West, M. (2004). Sparse graphical models for exploring gene expression data. Journal of Multivariate Analysis, 90(1):196–212. Dominici, F., Zeger, S. L., and Samet, J. M. (2000). A measurement error model for time-series studies of air pollution and mortality. Biostatistics, 1(2):157–175. Eichler, M. (2012). Causality, chapter Causal inference in time series analysis, 327–354. Wiley Series in Probability and Statistics. Fourer, R., Gay, D., and Kernighan, B. W. (2002). The AMPL book. Duxbury Press, Pacific Grove. George, E. I. (2000). The variable selection problem. Journal of the American Statistical Association, 95(452):1304–1308. George, E. I. and McCulloch, R. E. (1997). Approaches for Bayesian variable selection. Statistica sinica, 7(2):339–373. Gorrostieta, C., Fiecas, M., Ombao, H., Burke, E., and Cramer, S. (2013). Hierarchical vector auto-regressive models and their applications to multi-subject effective connectivity. 18 Griffin, M. P., Lake, D. E., and Moorman, J. R. (2005). Heart rate characteristics and laboratory tests in neonatal sepsis. Pediatrics, 115(4):937–941. Hastie, T. and Efron, B. (2013). Least Angle Regression, Lasso and Forward Stagewise. http://cran.r-project.org/web/packages/lars/lars.pdf. Haufe, S., Nolte, G., Mueller, K.-R., and Kr¨amer, N. (2010). Sparse causal discovery in multivariate time series. JMLR W&CP, 6:97–106. Hsu, N.-J., Hung, H.-L., and Chang, Y.-M. (2008). Subset selection for vector autoregressive processes using Lasso. Computational Statistics & Data Analysis, 52(7):3645–3657. Hu, J. and Hu, F. (2009). Estimating equation–based causality analysis with application to microarray time series data. Biostatistics, 10(3):468–480. Kalli, M. and Griffin, J. E. (2014). Time-varying sparsity in dynamic regression models. Journal of Econometrics, 178(2):779–793. Kojima, K., Yamaguchi, R., Imoto, S., Yamauchi, M., Nagasaki, M., Yoshida, R., Shimamura, T., Ueno, K., Higuchi, T., Gotoh, N., et al. (2009). A state space representation of VAR models with sparse learning for dynamic gene networks. Genome Informatics, 22:56–68. Lee, J. and Leyffer, S. (2012). Mixed integer nonlinear programming. Springer. Li, J. (2012). Monetary policy analysis based on Lasso-Assisted Vector Autoregression (Lavar). Available at SSRN 2017877. Lozano, A. C., Abe, N., Liu, Y., and Rosset, S. (2009). Grouped graphical Granger modeling for gene expression regulatory networks discovery. Bioinformatics, 25(12):i110–i118. L¨utkepohl, H. (2005). New introduction to multiple time series analysis. Springer Science & Business Media. Shojaie, A. and Michailidis, G. (2010). Discovering graphical granger causality using the truncating Lasso penalty. Bioinformatics, 26(18):i517–i523. Song, S. and Bickel, P. J. (2011). Large vector auto regressions. arXiv preprint arXiv:1106.3915. Tibshirani, R. (1996). Regression shrinkage and selection via the Lasso. Journal of the Royal Statistical Society. Series B (Methodological), 267–288. 19 Vald´es-Sosa, P. A., S´anchez-Bornot, J. M., Lage-Castellanos, A., Vega-Hern´andez, M., Bosch-Bayard, J., Melie-Garc´ıa, L., and Canales-Rodr´ıguez, E. (2005). Estimating brain functional connectivity with sparse multivariate autoregression. Philosophical Transactions of the Royal Society B: Biological Sciences, 360(1457):969–981. Yang, Y. and Zou, H. (2015). Group Lasso Penalized Learning Using A Unified BMD Algorithm. https://cran.r-project.org/web/packages/gglasso/gglasso.pdf. Yuan, M. and Lin, Y. (2006). Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(1):49–67. Zhao, P., Rocha, G., and Yu, B. (2009). The composite absolute penalties family for grouped and hierarchical variable selection. The Annals of Statistics, 3468– 3497. Zou, H. (2006). The adaptive Lasso and its oracle properties. Journal of the American Statistical Association, 101(476):1418–1429. Abbreviation Method Expression VAR Classic vector autoregressive approach (1) Lasso-VAR Vector autoregression + Lasso (2) Group Lasso Vector autoregression + Group Lasso + 5-fold cross-validation (3) SC-VAR Sparsity-controlled vector autoregressive method (P) Table 1: Summary of methods under comparison in the computational illustrations Abbreviation Name N Ttrain Ttest pReference Google flu Google Flu Trends 48 221 220 1 Davis et al. (2012) Air pollution Concentration levels of air pollutants 5 4185 4185 4 Davis et al. (2012) Table 2: Summary of real databases used for numerical illustrations 20 Real graph 1 23 4 5 0.7 0.4 0.5 0.4 0.4 0.4 1 23 4 5 0.5 1 23 4 5 0.5 1 23 4 5 Naive solution , MSE =2.36 1 23 4 5 0.7 1 23 4 5 1 23 4 5 1 23 4 5 Lasso solution, MSE =2.63 1 23 4 5 0.12 1 23 4 5 1 23 4 5 1 23 4 5 SC−VAR solution , MSE =1.99 1 23 4 5 1 23 4 5 0.81 1 23 4 5 1 23 4 5 Figure 1: Graphical representation of a simulated VAR and its sparse counterparts. Naive algorithm (center) versus the proposed SC-VAR method (right). 21 Density 10% Density 50% Density 90% Lasso SC-VAR Lasso SC-VAR Lasso SC-VAR VT= 1 1.09 1.04 1.23 1.15 1.23 1.17 VT= 2 1.03 1.00 1.15 1.07 1.20 1.11 VT= 3 1.01 0.99 1.09 1.03 1.15 1.08 VT= 4 0.99 0.99 1.05 1.02 1.12 1.05 VT= 5 0.99 0.99 1.03 1.01 1.09 1.04 VT= 6 1.00 0.99 1.01 1.01 1.06 1.04 VT= 7 1.00 0.99 1.00 1.01 1.03 1.03 VT= 8 1.00 0.99 1.00 1.01 1.02 1.03 VT= 9 1.00 0.99 1.00 1.01 1.00 1.03 VT= 10 1.00 0.99 1.00 1.01 1.00 1.03 Table 3: Normalized median MSE under the Lasso and the SC-VAR for VAR(1) processes generated with a 10, 50 and 90% density, for different requirements of sparsity (VT=1,...,10) 22 Density 10% 50% 90% Real VAR MSE =1 MSE =1 MSE =1 VA= 10 SC-VAR MSE =1.19 MSE =1.32 MSE =1.78 VT= 1 Lasso MSE =1.54 MSE =1.9 MSE =2.4 SC-VAR MSE =1.34 MSE =1.35 MSE =1.8 VT= 10 Lasso MSE =1 MSE =1 MSE =1 SC-VAR MSE =0.99 MSE =1 MSE =1.03 Figure 2: Real graphs of randomly selected instances of processes with different densities, represented along with their VAR solutions. SC-VAR and Lasso solutions are represented for a couple of levels of sparsity, VT= 1 and VT= 10, and the SC-VAR is also depicted fixing solely V A=10. 23 Density 10% Density 50% Density 90% Group Lasso SC-VAR Group Lasso SC-VAR Group Lasso SC-VAR VS= 1 1.13 0.95 1.11 1.01 1.12 1.02 VS= 2 1.15 0.90 1.11 0.96 1.13 0.99 VS= 3 1.16 0.90 1.14 0.90 1.14 0.94 VS= 4 1.17 0.90 1.17 0.89 1.15 0.91 VS= 5 1.18 0.90 1.20 0.87 1.18 0.89 VS= 6 1.19 0.90 1.24 0.87 1.21 0.88 VS= 7 1.20 0.90 1.28 0.87 1.25 0.87 VS= 8 1.20 0.90 1.30 0.87 1.31 0.87 VS= 9 1.21 0.90 1.33 0.87 1.40 0.87 VS= 10 1.24 0.90 1.36 0.87 1.48 0.87 Table 4: Normalized median MSE under the Group Lasso and the SC-VAR for VAR(2) processes generated with a 10, 50 and 90% density, for different requirements of sparsity in terms of number of causal features (VS=1,...,10) 24 ●● ● ●●●●●● ● 2 4 6 8 10 0.5 1.0 1.5 2.0 VS Elapsed time in seconds ●● ● ●●●●●● ● 2 4 6 8 10 0.5 1.0 1.5 2.0 VS Elapsed time in seconds 2 4 6 8 10 0.5 1.0 1.5 2.0 VS Elapsed time in seconds 2 4 6 8 10 0.5 1.0 1.5 2.0 VS Elapsed time in seconds 2 4 6 8 10 0.5 1.0 1.5 2.0 VS Elapsed time in seconds 2 4 6 8 10 0.5 1.0 1.5 2.0 VS Elapsed time in seconds ●density =0.1 density =0.5 density =0.9 p=1 p=2 Figure 3: Median elapsed times in seconds taken by the SC-VAR to be solved for different requirements of the sparsity in VAR(1) and VAR(2) processes with various levels of density. 0 10 20 30 40 0.5 1.0 1.5 2.0 2.5 3.0 Google flu database VT MSE 0 10 20 30 40 0.5 1.0 1.5 2.0 2.5 3.0 SC−VAR Lasso Figure 4: Normalized MSEs for the SC-VAR and the Lasso under different levels of sparsity. 25