LASSO-Driven Inference in Time and Space

LASSO-Driven Inference in Time and Space

Introduction

Many applications in statistics, economics, finance, biology and psychology are concerned with a system of ultra high-dimensional objects that communicate within complex dependency channels. Given a complex system involving many factors, one builds a network model by taking a large set of regressions, i.e. regressing every factor in the system on a large subset of other factors. Examples include analysis of financial systemic risk by quantile predictive graphical models with LASSO hautsch2014financial,hardle2016tenet,BCC2016, limit order book network modeling via the penalized vector autoregressive approach chen_2018, analysis of psychology data with temporal and cross- sectional dependencies epskamp2018gaussian. Another example is quantifying the spillover effects or externalities for a social network, especially when the social interactions (or the interconnectedness) is not obvious manresa2013estimating. Besides, there are numerous applications concerning association network analysis in other fields of applied statistics; see Chapter 7 in kolaczyk2014statistical. In general, a step-by-step LASSO procedure is very helpful for the correlation network formation. In pursuing a highly structural approach, one certainly favors a simple set of regressions that allows multiple insights on the statistical structure of the data. Therefore, a sequence of regressions with LASSO is a natural path to take. Especially in cases of reduced forms of simultaneous equation models and structural vector autoregressive models, one can attain valuable pre-information on the core structure by running a set of simple regressions with LASSO shrinkage.

A first important question arising in this framework is how to decide on a unified level of penalty. In this article we advocate an approach to selecting the overall level of the tuning parameter in a system of equations after performing a set of single step regressions with shrinkage. A feasible (block) bootstrap procedure is developed and the consistency of parameter estimation is studied. In addition, we provide a uniform near-oracle bound for the joint estimators. The proposed technique is applicable to ultra-high dimensional systems of regression equations with high-dimensional regressors.

A second crucial issue is to establish simultaneous inference on parameters, which is an important question regarding network topology inference. For example, in a large-scale linear factor pricing model, it is of great interest to check the significance of the intercepts of cross sectional regressions (connected with zero pricing errors), e.g. pesaran2017testing. Our approach is an alternative testing solution compared to the Wald test statistics proposed therein. To achieve the goal of simultaneous inference, we develop a uniform robust post-selection or post-regularization inference procedure for time series data. This method is generated from a uniform Bahadur representation of de-biased instrumental variable estimators. In particular, we need to establish maximal inequalities for empirical processes for a general Huber's $Z$-estimation. Note that the commonly used technique for independent data, such as the symmetrization technique, is not directly applicable in the dependent data case; see Chapter 11.6 of kosorok2008introduction for a related overview.

Our contribution lies in three aspects. First, we select the penalty level by controlling the aggregated errors in a system of high-dimensional sparse regressions, and we establish the bounds on the estimated coefficients. Furthermore, we show the implication of the restricted eigenvalue (RE) condition at a population level. Secondly, an easily implemented algorithm for effective estimation and inference is proposed. In fact, the offered estimation scheme allows us to make local and global inference on any set of parameters of interest. Thirdly, we run numerical experiments to illustrate good performance of our joint penalty relative to the single equation estimation, and we show the finite sample improvement of our multiplier block bootstrap procedure on the parameter inference. Finally, an application of textual sentiment spillover effects on the stock returns in a financial market is presented.

In the literature, the fundamental results on achieving near oracle rate for penalized $\ell_1$-norm estimators are developed by bickel2009simultaneous. There are many related articles on deriving near-oracle bounds using the $\ell_1$-norm penalization function for the i.i.d. case, such as BCH2011,belloni2009least. There are also many extensions to the LASSO estimation with dependent data. For example, {basu2015regularized study the consistency of the estimator in sparse high-dimensional Gaussian time series models; kock2015oracle consider the high-dimensional near-oracle inequalities in large vector autoregressive (VAR) models; lin2017regularized look at the regularized estimation and testing for high-dimensional multi-block VAR models.} However, the majority of the literature imposes a {Gaussian} or sub-Gaussian assumption on the error distribution; this is rather restrictive and excludes heavy tail distributions. For dependent data, wu2016performance discuss the possibility of relaxing the sub-Gaussian assumption by generalizing Nagaev-type inequalities allowing for only moment assumptions. For the case of LASSO the analysis assumes the fixed design, which rules out the most important applications mentioned earlier in the introduction.

Theoretically, the LASSO tuning parameter selection requires characterizing the asymptotic distribution of the maximum of a high dimensional random vector. CCK13AoS develop a Gaussian approximation for the maximum of a sum of high-dimensional random vectors, which is in fact the basic tool for modern high-dimensional estimation. Here it is applied to the LASSO inference. Moreover, CCK13testing deliver results for the case of $\beta$-mixing processes. Although it is quite common to assume a mixing condition {which is at base a concept yielding asymptotic independence}, it is not in general easy to verify the condition for a particular process, and some simple linear processes can be excluded from the strong mixing class, andrews1984non. With an easily accessible dependency concept, ZW15gaussian derive Gaussian approximation results for a wide class of stationary processes. Note that the dependence measure is linked to martingale decompositions and is therefore readily connected with a pool of results on tail probabilities, moment inequalities and central limit theorems of martingale theory. Our results are built on the above-mentioned theoretical works and we extend them substantially to fit into the estimation in a system of regression equations. In particular, our LASSO estimation is with random design for dependent data; therefore, we need to deal with the population implications of the Restricted Eigenvalue (RE) condition. {Moreover, we show the interaction between the tail assumption and the dimensionality of the covariates in our theoretical results.}

In the meantime, the issue of simultaneous inference is challenging and has motivated a series of research articles. For the case of i.i.d. data, BCH2011,BCH2014, zhang2014debiased, JM2014, van2014desparsified, neykov2015unified, DML, zhu2017linear, among others, develop confidence intervals of low-dimensional variables in high-dimensional models with various forms of de-biased/orthogonalization methods. Still in the case of i.i.d. data, BCK15Bio establish a uniform post-selection inference for the target parameters defined via de-biased Huber's $Z$-estimators when the dimension of the parameters of interest is potentially larger than the sample size, where they employ the multiplier bootstrap to the estimated residuals. Wild and residual bootstrap-assisted approaches are also studied in dezeure2016high,zhang2017simultaneous for the case of mean regression. {And more recently, krampe2018 extend the approaches to test large groups of coefficients in sparse VAR models.} We pick up the line of the inference analysis of BCK15Bio and employ it in a temporal and cross-sectional dependence framework, thus making it applicable to a rich class of high-dimensional time series. {This allows us to embed the high-dimensional VAR model as a special case.} Our core proof strategy is different, as it is well known that the technique for handling the suprema of empirical processes indexed by functional classes with dependent data is not the same as in i.i.d. cases. For instance, the key Bahadur representation in BCK15Bio applies maximal inequalities derived in chernozhukov2014gaussian for i.i.d. random variables, while we derive the key concentration inequalities based on a martingale approximation method.

{Our proposed estimation framework is complement to the literature on model selection for Gaussian Graphical model (GGM) (see e.g. yuan2007model), which has a wide spectrum of applications in statistics. A GGM can be connected with LASSO regression for estimating sparse correlation networks, and therefore is equivalent to our context with a partial correlation network, meinshausen2006high. In particular, we may find an equation-by-equation relationship to the GGM, and we acknowledge that a similar framework with spatial temporal dependence can be developed. In addition, there is a big literature on social network analysis, which embeds the network information into a dynamic model in advance; see for example zhu2017network,zhu2016network,chen2019tail,huang2016statistical. Relatively, our approach is less structural as we treat the network structure to be unknown and uncover it using LASSO. }

The following notations are adopted throughout this paper. For a vector $v=(v_1,\ldots,v_p)^\top$, let $|v|_\infty\stackrel{\mathrm{def}}{=}\max_{1\leqslant j\leqslant p}|v_j|$ and $|v|_s\stackrel{\mathrm{def}}{=}(\sum_{j=1}^p|v_j|^s)^{1/s}$, $s\geqslant1$. For a random variable $X$, let $\|X\|_q\stackrel{\mathrm{def}}{=}(\mathop{\mbox{\sf E}}|X|^q)^{1/q}$, $q>0$. For any function on a measurable space $g:\mathcal{W}\rightarrow\rm I\!R$, $\mathop{\mbox{\sf E}}_n(g)\stackrel{\mathrm{def}}{=} n^{-1}\sum_{t=1}^n\{g(\omega_t)\}$ and $G_n(g)\stackrel{\mathrm{def}}{=} n^{-1/2}\sum_{t=1}^n[g(\omega_t)-\mathop{\mbox{\sf E}}\{g(\omega_t)\}]$. Given two sequences of positive numbers $a_n$ and $b_n$, write $a_n\lesssim b_n$ if there exists constant $C>0$ (does not depend on $n$) such that $a_n/b_n\leqslant C$. For a sequence of random variables $x_n$, we use the notation $x_n\lesssim_{\operatorname{P}} b_n$ to denote $x_n=\mathcal{O}_{\operatorname{P}}(b_n)$. {For any {finitely discrete} measure $\mathcal Q$ on a measurable space, let $\mathcal L^q(\mathcal Q)$ denote the space of all measurable functions $f:Z\rightarrow{\rm I\!R}$ such that $\|f\|_{\mathcal Q,q}\stackrel{\mathrm{def}}{=}(\mathcal Q|f|^q)^{1/q}<\infty$, where $\mathcal Qf\stackrel{\mathrm{def}}{=}\int fd\mathcal Q$. For a class of measurable functions $\mathcal F$, the $\epsilon$-covering number with respect to the $\mathcal L^q(\mathcal Q)$-semimetric is denoted as $\mathcal N(\epsilon,\mathcal F, \|\cdot\|_{\mathcal Q,q})$, and let $\operatorname{ent}(\epsilon,\mathcal F)=\log\sup_{\mathcal Q}\mathcal N(\epsilon\|\bar F\|_{\mathcal Q,q},\mathcal F,\|\cdot\|_{\mathcal Q,q})$ with $\bar F=\sup_{f\in\mathcal F}|f|$ (the envelope) denote the uniform entropy number. It should be noted that we suppress the notation of the outer expectation $\mathop{\mbox{\sf E}}^*$ to $\mathop{\mbox{\sf E}}$ and outer probability $\operatorname{P}^*$ to $\operatorname{P}$ when measurability issues are encountered. Details may be found in the Chapter 1 of van1996weak.}

The rest of the article is organized as follows. Section (ref) shows the system model with a few examples. Section (ref) introduces the sparsity method for effective prediction and provides an algorithm for the joint penalty level of LASSO via bootstrap. In Section (ref) we propose approaches to implementing individual and simultaneous inference on the coefficients. Main theorems are listed in Section (ref). In Section (ref) and (ref) we deliver the simulation studies and an empirical application on textual sentiment spillover effects. The technical proofs and other details are given in the supplementary materials. The codes to implement the algorithms are publicly accessible via the website \href{https://github.com/QuantLet/LASSO_Time_Space}{\raisebox{-1pt}\, www.quantlet.de}.

The System Model

In this section, we present a general framework which covers many applications in statistics. Consider the system of regression equations (SRE): $$ Y_{j,t} = X_{j,t}^\top\beta_{j}^0 + {\varepsilon}_{j,t}, \quad \mathop{\mbox{\sf E}} {\varepsilon}_{j,t} X_{j,t} = 0, \quad j=1,...,J, \quad t=1,\ldots, n, $$ where $X_{j,t}=(X_{jk,t})_{k=1}^{K_j}$. {Without loss of generality, we assume the dimension of the covariates is identical among all equations thereafter, namely $K_j = \dim(X_{j,t})\equiv K$, for $j=1,\ldots,J$.} We allow the dimension $K$ of $X_{j,t}$ and the number of equations, $J$ to be large, potentially larger than $n$, which creates an interplay with the tail assumptions on the error processes ${\varepsilon}_{j,t}$. Both spatial and temporal dependency are allowed and we will obtain results on prediction and inference.

The SRE framework is a system of regression equations, which includes the following important special cases.

{We refer to Section (ref) in the supplementary materials for more practical examples.}

Effective Prediction Using Sparsity Method

In this section, we present our model setup and the LASSO estimation algorithm, including the joint penalty selection procedure.

Sparsity in SRE

The general SRE structure makes it possible to predict $Y_{j,t}$ using $X_{j,t}$ effectively. Note that the dimension of $X_{j,t}$ is large, potentially larger than $n$. Without loss of generality we assume exact sparsity of $\beta_{j}^0$ throughout the paper:

where the $\ell_0$-norm, $|\cdot|_0$, is the number of nonzero components of a vector.

For this situation one employs an $\ell_1$-penalized estimator of $\beta_{j}^0$ of the form:

where $\lambda$ is the joint "optimal" penalty level and $\Psi_{jk}$'s are penalty loadings, which are defined below in ((ref)).

A first aim is to obtain performance bounds with respect to the prediction norm: $$ |\widehat \beta_j - \beta_{j}^0|_{j, pr} \stackrel{\mathrm{def}}{=} \bigg[\frac{1}{n} \sum_{t=1}^n \big\{ X_{j,t}^\top(\widehat \beta_j - \beta_{j}^0)\big\}^2 \bigg]^{1/2}, $$ {where the outside $j$ indicates to use the covariates in the $j$th equation $X_{j,t}$ in computing the prediction norm,} and the Euclidean norm: $$ |\widehat \beta_j - \beta_{j}^0|_{2} \stackrel{\mathrm{def}}{=} \big\{\sum_{k=1}^K(\widehat \beta_{jk} - \beta_{jk}^0)^2 \big\}^{1/2}. $$ To achieve good performance bounds, we first consider "ideal" choices (IC) of the penalty level and the penalty loadings. Let $$ S_{jk} = \frac{1}{\sqrt{n}} \sum_{t=1}^n {\varepsilon}_{j,t} X_{jk,t}, $$ where for a moment we assume to be able to observe ${\varepsilon}_{j,t} = Y_{j,t} - X_{j,t}^\top\beta_{j}^0$. In practice one obtains an approximation by stepwise LASSO. Set

where $c >1$, e.g., $c=1.1$, and $1-\alpha$ is a confidence level, e.g. $\alpha =0.1$, {where the long run variance is denoted by $\operatorname{avar}$}.

Theoretically, we can characterize the rate of $\lambda^0(1-\alpha)$ by the tail probability of $S_{jk}$ (see Theorem (ref)), also via Gaussian Approximation as in corollary (ref). To calculate $\lambda^0(1-\alpha)$ from data, we can also use a Gaussian approximation based on: $$ Q(1-\alpha) \stackrel{\mathrm{def}}{=} (1-\alpha)-\textrm{quantile of } {2c \sqrt{n}}\max_{1 \leqslant j\leqslant J, 1\leqslant k \leqslant K}|Z_{jk}/\Psi_{jk}|, $$ where $\{Z_{jk}\}$ are multivariate Gaussian centered random variables with the same {long run covariance} structure as {$\{S_{jk}\}$}. Alternatively, we can employ a multiplier bootstrap procedure to estimate IC empirically to achieve a better finite sample performance; see for example CCK13AoS. In case of dependent observations over time, it is understood that data cannot be resampled directly as in the the i.i.d. case, as the dependency structure of the underlying processes will be lost. A usual solution to this problem is to consider a block bootstrap procedure, where the data are grouped into blocks, resampled and concatenated. In particular, we will adopt an estimate of IC by a multiplier block bootstrap procedure. {The theoretical properties of LASSO and the tuning parameter choices are presented in Section (ref)-(ref).}

Multiplier Bootstrap for the Joint Penalty Level

In this subsection, we introduce an algorithm to approximate the joint penalty level via a block multiplier bootstrap procedure, which is particularly non-overlapping block bootstrap (NBB). Consider the system of equations with dependent data:

The bootstrap consistency regarding $Z^{[B]}_{jk}$ is proved in Theorem (ref).

Valid Inference on the Coefficients

With a reasonable fitting of LASSO on hand, we can proceed to investigate the issue of simultaneous inference. This section focuses on SRE of Example 2. We allow the covariates in each equation to be different.

The basic idea to facilitate inference is to formulate the estimation in a semi-parametric framework. With partialing out the effect of the nonparametric coefficient(s), we can achieve the desired estimation accuracy of the parametric component of interest. This trick is referred to as "Neyman orthogonalization". Notably, the procedure is equivalent to the well known de-sparsification procedure in the mean square loss case, which is developed for the inference on the estimated zero coefficients by LASSO. It thus serves the same purpose of generating a (robust) de-sparsified estimation for LASSO inference.

We list three algorithms to estimate $\beta^0_{jk}$. Algorithm 1 is easy to implement and algorithm 2 is tailored to the cases of heavy-tailed distribution of the error term, as Least Absolute Deviation (LAD) regression is well known to be robust against outliers. Algorithm 3 considers a double selection procedure aimed at remedying the bias due to omitted variables by one step selection, while also accounting for the cases of heteroscedastic errors.

Algorithm 1: LS-based algorithm

Algorithm 2: LAD-based algorithm

{The theoretical properties of the estimators $\widehat\beta_{j(-k)}^{[1]}$ and $\widehat\gamma_{j(-k)}$ in S1 and S2 are provided in Corollary (ref) or (ref) (see Corollary (ref) or (ref) in the supplementary correspondingly if the joint penalty over equations is employed), and Theorem (ref) for post LASSO, respectively. The uniform Bahadur representation and the Central Limit Theorem of the estimator $\widehat\beta^{[2]}_{jk}$ in S3 or S3$'$ are established in Theorem (ref) and Corollary (ref).}

Algorithm 3: Double selection-based algorithm

As shown in BCH2011 and BCK15_sup, the double selection approach in S3$''$ or S3$'''$ creates an orthogonality condition with respect to the space spanned by the covariates selected by both steps, and thus generates an orthogonal relation to any space spanned by a linear projection of the covariates, e.g. $\widehat v_{jk,t}$. Therefore, the inference on the parameters may still be applied as in the framework of Algorithm \hyperref[algo1]{1} and \hyperref[algo2]{2}. {Therefore, one may still find the theoretical properties of estimators in S1$''$, S2$''$, S3$''$ (S3$'''$) in Section (ref) according to the links mentioned above.}

Confidence Interval for a Single Coefficient

We discuss an inference framework developed for a single coefficient obtained from the aforementioned algorithms.

Let $\psi_{jk}(Z_{j,t},\beta_{jk},h_{jk})$ denote the score function, where $Z_{j,t}=(Y_{j,t},X^\top_{j,t})^\top$, $h_{jk}(X_{j(-k),t})=(X_{j(-k),t}^\top\beta_{j(-k)},X_{j(-k),t}^\top\gamma_{j(-k)})^\top$. Consider the LAD-based case with $\psi_{jk}(Z_{j,t},\beta_{jk},h_{jk})=\{1/2-\boldsymbol{1}(Y_{j,t}\leqslant X_{jk,t}\beta_{jk}+X_{j(-k),t}^\top\beta_{j(-k)})\}v_{jk,t}$, define $\omega_{jk} \stackrel{\mathrm{def}}{=}\mathop{\mbox{\sf E}} \{(\frac{1}{\sqrt{n}}\sum_{t=1}^n\psi^0_{jk,t})^2\}=\sum_{\ell=-(n-1)}^{n-1}(1-\frac{|\ell|}{n})\operatorname{cov}(\psi^0_{jk,t},\psi^0_{jk,(t-\ell)})$ with $\psi_{jk,t}^0\stackrel{\mathrm{def}}{=}\psi_{jk}(Z_{j,t},\beta^0_{jk},h^0_{jk})$, and $\phi_{jk}\stackrel{\mathrm{def}}{=}\frac{\partial \mathop{\mbox{\sf E}}\{\psi_{jk}(Z_{j,t},\beta,h^0_{jk})\}}{\partial \beta}|_{\beta=\beta_{jk}^0}$.

Suppose we are interested in testing $H_0: \beta^0_{jk}=0$. For this purpose we employ the uniform Bahadur representation (Theorem (ref)) to construct the confidence interval via a multiplier bootstrap procedure. In particular, the distribution of the asymptotically pivotal statistics:

is approximated via its block multiplier bootstrap counterpart:

where {$\widehat\zeta_{jk,t}$ are pre-estimators of $\zeta_{jk,t}=-\phi^{-1}_{jk}\sigma_{jk}^{-1}\psi^0_{jk,t}$ such that \\ $\underset{(j,k),(j',k')}{\max}|\sum_{i=1}^{l_n} \widehat\eta_{j'k',i}\widehat\eta_{jk,i} - \sum_{i=1}^{l_n} \eta_{j'k',i}\eta_{jk,i}|=\mbox{\tiny $\mathcal{O}$}_{\operatorname{P}}(\{\log(JK)\}^{-2})$, with $\eta_{jk,i} \stackrel{\mathrm{def}}{=} \frac{1}{\sqrt{n}} \sum_{l=(i-1)b_n+1}^{ib_n} \zeta_{jk,l} $ and $\widehat{\eta}_{jk,i} \stackrel{\mathrm{def}}{=}\frac{1}{\sqrt{n}} \sum_{l=(i-1)b_n+1}^{ib_n} \widehat \zeta_{jk,l}$}, $e_{j,i}$ are independently drawn from $\operatorname{N}(0,1)$, $l_n$ and $b_n$ are the numbers of blocks and block size, respectively. {More discussion on how one can construct the consistent pre-estimators $\widehat\zeta_{jk,t}$ is stated in the supplementary material; see Comment (ref).}

Let $\widehat{\sigma}_{jk}$ be any consistent estimator of $\sigma_{jk}$. Then the confidence interval is given by

where $q^{\ast}_{jk}(1-\alpha)$ is the $(1-\alpha)$ quantile of the bootstrapped distribution of $|T_{jk}^{\ast}|$.

{

}

Joint Confidence Region for Simultaneous Inference

We now continue to extend the single coefficient inference to simultaneous inference on a set of coefficients. As shown in the practical examples in Section (ref), it is essential to conduct simultaneous inference on a group of parameters $G$. In this case, the null hypothesis is: $\mathbf{H}_0: \beta_{jk}^0 = 0$, $\forall (j,k) \in G $, and the alternative $\mathbf{H}_A: \beta_{jk}^0 \neq 0$, for some $(j,k) \in G$, where the group $G$ is a set of coefficients with cardinality $|G|$. Suppose for the $j$-th equation there are $p_j$ target coefficients and the cardinality $|G| = \sum^J_{j=1} p_{j}$. This can be understood as a multiple estimation problem compared to Section (ref). Without loss of generality, we can rearrange the order of the variables and rewrite the regression equation for each $j$ as (consider the LAD-based model here)

One follows the algorithms to obtain $\widehat{\beta}_{jl} (1\leqslant l\leqslant p_{j})$ for each $j$. Then the idea of simultaneous inference is very straightforward. We aggregate the statistics $T_{jk}$ in (ref) by taking the maximum and minimum over the set $G$. Finally, the component-wise confidence interval is constructed with the quantiles of the bootstrap statistics over all bootstrap samples.

Denote $q^{\ast}_{G}(1-\alpha)$ as the $(1-\alpha)$ quantile of $\underset{(j,k)\in G}{\max} |T^{\ast}_{jk}|$. A joint confidence region is then:

and for each component $(j,k)\in G$, the confidence interval $\widetilde{\operatorname{CI}}^\ast_{jk}(\alpha)$ is given by $ [\widehat{\beta}^{[2]}_{jk}-\widehat{\sigma}_{jk} n^{-1/2}q^{\ast}_{G}(1-\alpha), \widehat{\beta}^{[2]}_{jk}+\widehat{\sigma}_{jk}n^{-1/2}q^{\ast}_{G}(1-\alpha)]$. We show in Corollary (ref) the consistency of this bootstrap confidence band for simultaneous inference. {Note that when there is only one parameter in $G$ for inference, the joint confidence region (ref) will reduce to the single parameter confidence interval (ref) as a special case.}

Main Theorems

In this section, we present the theoretical foundations for the procedures given earlier. In particular, we discuss the properties of the theoretical choices of penalty level and the validity of the other two empirical choices, as well as the theoretical support for the simultaneous inference.

Throughout the whole section, we define $S_{jk} \stackrel{\mathrm{def}}{=} n^{-1/2} \sum_{t=1}^n \varepsilon_{j,t} X_{jk,t}$, $S_{j\cdot}= (S_{jk})_{k=1}^K$, and $\Psi_{jk}\stackrel{\mathrm{def}}{=}\sqrt{\operatorname{avar}(S_{jk})}$, which is the square root of the long-run variance of $X_{jk,t}{\varepsilon}_{j,t}$, namely \\$\{\sum^{\infty}_{\ell=-\infty} \mathop{\sf E} (X_{jk,t}X_{jk,(t-\ell)}{\varepsilon}_{j,t} {\varepsilon}_{j,(t-\ell)} )\}^{1/2}$. Recall that for a single equation LASSO, we select the penalty in the following ways:

Near Oracle Inequalities under IC

We first provide the near oracle inequalities for the single equation LASSO estimation $\tilde{\beta}_j$ obtained from (ref) under the ideal choices (IC). For this purpose, a few assumptions and definitions are required.

It should be noted that \hyperref[A1]{(A1)} admits a wide class of processes. The largest value of $\varsigma$ which ensures a finite dependence adjusted norm characterizes the dependency structure of the process. The moment-based measure is directly connected with the impulse functions. {A few examples for univariate time series $Z_t$ are listed in Appendix (ref) in the supplementary materials.}

{Lemma (ref) follows Theorem 1 of belloni2009least. As the proof is built on inequalities and for the case of dependent data \hyperref[A1]{(A1)} they remain unchanged, we omit the detailed proof here.} To further characterize the rate of IC, we provide a tail probability for $2 c\sqrt{n} \underset{1\leqslant k \leqslant K}{\max}|S_{jk}{/\Psi_{jk}}|$ under the moment assumption \hyperref[A3]{(A3)}. In particular, the rate depends on the dependence adjusted norm $\|X_{jk,\cdot}{\varepsilon}_{j,\cdot}\|_{q,\varsigma}$.

Under the choice (IC) $\lambda_j^0(1-\alpha)$ is given by the $(1-\alpha)$ quantile of $2 c\sqrt{n} \underset{1\leqslant k \leqslant K}{\max}|S_{jk}{/\Psi_{jk}}|$, combining the results of Lemma (ref) and Theorem (ref) we can get the bounds for $\lambda_j^0(1-\alpha)$ and further obtain the oracle inequalities as in Corollary (ref).

We now discuss the case of sub-Gaussian tail or sub-exponential tail, which is mostly assumed in the literature.

Gaussian Approximation for Dependent Data

Now we look at the validity of the choice of $Q_j(1-\alpha)$, which relies on a Gaussian approximation theorem. First we define the Kolmogorov distance between any two $K$-dim random vectors.

{For each single equation $j$, aggregate the dependence adjusted norm over $k=1,\ldots,K$:

where $q\geqslant1$ and $\varsigma>0$. Moreover, define the following quantities

} {It is worth noting that the norm $ \||X_{j,\cdot}|_\infty\|_{q,\varsigma}$ is a kind of aggregated dependence adjusted norm for a vector of processes in comparison to the dependence adjusted norm for a univariate process as in Definition (ref). }

Some additional assumptions are required. Define $L_{1,j} = \{\Phi_{j,4,\varsigma}\Phi_{j,4,0} (\log K)^2\}^{1/\varsigma}$, $W_{1,j}$ $=$ $(\Phi^6_{j,6,0}+ \Phi^4_{j,8,0})\{\log(Kn)\}^7$, $W_{2,j}$ $=$ $\Phi^2_{j,4,\varsigma}\{\log(Kn)\}^4$, $W_{3,j}$ $=$ $[n^{-\varsigma} \{\log (Kn)\}^{3/2} \Theta_{j,2q,\varsigma}]^{1/(1/2-\varsigma-1/q)}$, $N_{1,j} =(n/\log K)^{q/2} \Theta_{j,2q, \varsigma}^{q}$, $N_{2,j}=n(\log K)^{-2}\Phi_{j,4,\varsigma}^{-2}$, $N_{3,j} = \{n^{1/2}(\log K)^{-1/2}\Theta^{-1}_{j,2q, \varsigma}\}^{1/(1/2-\varsigma)}$.

The assumptions impose mild restrictions on the dependency structure of covariates and error terms. They include a wide class of potential correlation and heterogeneity (including conditional heteroscedasticity), with possible allowance of the lagged dependent variables. {Two examples of large VAR and ARCH for high-dimensional time series can be found in Appendix (ref) in the supplementary materials.}

{

}

{

}

Theorem (ref) justifies the choice of $\lambda_j$ and $\tilde{Q}_j(1-\alpha)$, which leads to the following corollary:

{ It is worth noting that in practice the variance involved in the Gaussian approximation in (ref) is not known; we shall discuss how we estimate the variance and also the validity of the Gaussian approximation result with an estimated variance. Given the realization $X_{j,1}{\varepsilon}_{j,1},\ldots,X_{j,n}{\varepsilon}_{j,n}$, we propose to estimate the $K\times K$ long-run variance-covariance matrix $\Sigma_j$ for $j=1,\ldots,J$ as follows, given $\mathop{\mbox{\sf E}} X_{j,t}{\varepsilon}_{j,t}=0$, and consider:

Moreover, the following corollary ensures that the Gaussian approximation results still hold if we use the estimate in (ref).

}

{ It should be noted that given the Gaussian approximation results in Theorem (ref), we can have a refined bound for $\lambda_j^0(1- \alpha)$ and also the oracle inequalities under IC.

We note that the allowed dimension $K$ is still of polynomial rate restricted by \hyperref[A4]{(A4)}.

}

Multiplier Block Bootstrap Procedure

In this subsection, we discuss how $\Lambda_j(1-\alpha)$ is attainable via block bootstrap. The data over $t =1,\ldots,n $ are divided into $l_n$ blocks with the same number of observations $b_n$, $n=b_nl_n$ (without loss of generality), where $b_n, l_n\in\mathbb{Z}$.

Recall that $\Lambda_j(1-\alpha)=2c\sqrt{n}q^{[B]}_{j,(1-\alpha)}$, $q^{[B]}_{j,(1-\alpha)}$ is the $(1-\alpha)$ quantile of $\underset{ 1\leqslant k \leqslant K}{\max} |Z^{[B]}_{jk}{/\Psi_{jk}}|$, where $Z^{[B]}_{jk}$ are defined as

and $e_{j,i}$ are i.i.d. $\operatorname{N}(0,1)$ random variables independent of $X$ and ${\varepsilon}$.

In fact, the above construction relies on knowing the true residuals ${\varepsilon}_{j,t}$. In practice, one needs to pre-estimate them using a conservative choice of penalty levels and loadings. {We discuss the consistency rate of the bootstrap statistics with generated errors in the supplementary material; see Comment (ref) and Theorem (ref).}

{

Joint Penalty over Equations

Recall that the theoretical choice $\lambda^0(1-\alpha)$ is defined as the $(1-\alpha)$ quantile of $\underset{1\leqslant k\leqslant K, 1\leqslant j\leqslant J}{\max}2c\sqrt{n}|S_{jk}{/\Psi_{jk}}|$. The empirical choices of the joint penalty level can be:

Section (ref) in the supplementary material provides the main theorems for joint equation estimation. In particular, the dimension along $k=1,\ldots, K$ and $j=1,\ldots,J$ will be considered together by vectorization, resulting in the dimension of $KJ$. Following the results for the single equation (where $j$ is fixed), we generalize the theorems above to multiple equations case by changing the dimension from $K$ to $KJ$; see Section (ref) in the Appendix for more details.

}

Post-Model Selection Estimation

LASSO estimation is known to be biased especially for large coefficients. Therefore, a post-selection step helps to reduce the bias by running an OLS as a second step on the selected covariates in the first step. In particular, we consider the 2-step OLS post-LASSO estimator:

To provide the prediction performance bounds for the OLS post-LASSO estimators, we need the following restricted sparse eigenvalue (RSE) condition:

Here $p$ denotes the restriction on the length of the active set of $T_j^c$. When $T_j=\emptyset$, \hyperref[A5]{(A5)} is reduced to the standard sparse eigenvalue condition. Moreover, let $\mu_j(p)\stackrel{\mathrm{def}}{=}\sqrt{\phi_j(p)}/\tilde\kappa_j(p)$, and denote by $\widehat p_j\stackrel{\mathrm{def}}{=}|\widehat T_j\setminus T_j|$ the number of components outside $T_j\stackrel{\mathrm{def}}{=}\mathrm{supp}(\beta^0_j)=\{k\in\{1,\ldots,K\}: \beta_{jk}^0\neq0\}$ selected by LASSO in the first step.

{The performance bounds for the OLS post-LASSO estimator are shown in Theorem (ref) in the supplementary materials.}

Simultaneous Inference

This subsection develops theory corresponding to Section (ref). A key Bahadur representation which linearize the estimator for a proper application of the central limit theorem for inference is provided.

Recall that for each $j=1,\ldots,J$, the following model is considered

where we define $\gamma_{j(-k)}^0\stackrel{\mathrm{def}}{=}\arg\underset{\gamma_{j(-k)}}{\min}\mathop{\mbox{\sf E}}(X_{jk,t}- X_{j(-k),t}^{\top}\gamma_{j(-k)})^2$, {and let $F_{{\varepsilon}_j}$ denote the distribution function of ${\varepsilon}_{j,t}$.} In this subsection, we show the validity of the joint confidence region for simultaneous inference on $H_0: \beta_{jk}^0=0, \forall(j,k)\in G$, with $|G| = \sum^J_{j=1}p_j$. In particular, for $j=1,\ldots,J$, $\beta_{jk}^0\,(k=1,\ldots,p_j)$ are the target parameters. Theoretically, we formulate the estimation as a general $Z$-estimation problem, with the leading examples as the LAD/LS cases. Nevertheless, it can also include a more general class of loss functions.

For each $(j,k)\in G$, we define the score function as $\psi_{jk}\{Z_{j,t},\beta_{jk},h_{jk}(X_{j(-k),t})\}$, where $Z_{j,t}\stackrel{\mathrm{def}}{=}(Y_{j,t},X_{j,t}^\top)^\top$ and the vector-valued function $h_{jk}(\cdot)$ is a measurable map from ${\rm I\!R}^{K-1}$ to ${\rm I\!R}^M$ ($M$ is fixed). In particular, in our linear regression case we have $h_{jk}(X_{j(-k),t})=(X_{j(-k),t}^\top\beta_{j(-k)},X_{j(-k),t}^\top\gamma_{j(-k)})^\top$, and for the LAD regression $\psi_{jk}\{Z_{j,t},\beta_{jk},h_{jk}(X_{j(-k),t})\}=\{1/2-\boldsymbol{1}(Y_{j,t}\leqslant X_{jk,t}\beta_{jk}+X_{j(-k),t}^\top\beta_{j(-k)})\}(X_{jk,t}-X_{j(-k),t}^{\top}\gamma_{j(-k)})$.

Assume that there exists $s=s_n\geqslant1$ such that $|\beta_{j(-k)}^0|_0\leqslant s$, $|\gamma_{j(-k)}^0|_0\leqslant s$, for each $(j,k)\in G$. Moreover, we assume that the nuisance function $h^0_{jk}=(h^0_{jk,m})_{m=1}^M$ admits a sparse estimator $\widehat h_{jk}=(\widehat h_{jk,m})_{m=1}^M$ of the form $$ \widehat h_{jk,m}(X_{j(-k),t})=X_{j(-k),t}^\top\widehat\theta_{jk,m},\quad |\widehat\theta_{jk,m}|_0\leqslant s,\quad m=1,\ldots,M, $$ where the sparsity level $s$ is small compared to $n$ ($s\ll n$).

The true parameter $\beta_{jk}^0$ is identified as a unique solution to the moment condition

However, the object $\arg\,\underset{\beta_{jk} \in \widehat{\mathcal{B}}_{jk}}{\operatorname{zero}}\mathop{\mbox{\sf E}}_n|[\psi_{jk}\{Z_{j,t},\beta_{jk},h^0_{jk}(X_{j(-k),t})\}]|$ does not necessarily exist due to the discontinuity of the function $\psi_{jk}$. The estimator $\widehat\beta_{jk}$ is obtained as a $Z$-estimator by solving the sample analogue of (ref)

where $g_n \stackrel{\mathrm{def}}{=} \{\log(e|G|)\}^{1/2}$ and $\widehat{\mathcal{B}}_{jk}$ is defined in \hyperref[C2]{(C2)}.

We now lay out the following conditions needed in this section, which are assumed to hold uniformly over $(j,k)\in G$.

Conditions \hyperref[C1]{(C1)}-\hyperref[C4]{(C4)} and \hyperref[C7]{(C7)} assume mild restrictions on the $Z$-estimation problems. They include the LAD-based regression (used in Algorithm \hyperref[algo2]{2}) with non-smooth score function. {Conditions \hyperref[C2]{(C2)} and \hyperref[C8]{(C8)} imply that $\underset{(j,k)\in G}{\max}\|\psi^0_{jk,\cdot}\|_{2,\varsigma}\lesssim s^{1/2}\underset{f\in\mathcal F'}{\max}\|f(z_t)\|_{2}$ and $\big\|\underset{(j,k)\in G}{\max}|\psi^0_{jk,\cdot}|\big\|_{q,\varsigma}\lesssim s^{3/2} \|\bar F'(z_t)\|_q$.} In \hyperref[C5]{(C5)}, we suppose that the nuisance parameters have estimators with good sparsity and convergence rate properties. As discussed in previous sections, given the ideal choice of the tuning parameter, the oracle inequalities provided in Corollary (ref) {ensures that our proposed algorithms can produce the estimator of the form $| \widehat \beta_{j(-k)}^{[1]} - \beta_{j(-k)}^0|_{j, pr} \lesssim_{\operatorname{P}} \{\sqrt{s\log(a_n/\alpha)/n}\vee n^{1/q-1}(\varpi_na_n/\alpha)^{1/q}\}\underset{1\leqslant k \leqslant K}{\max}\|X_{jk,\cdot}{\varepsilon}_{j,\cdot}\|_{q,\varsigma}$, where for $\varsigma > 1/2-1/q$ (weak dependence case), $\varpi_n = 1$; for $\varsigma< 1/2-1/q$ (strong dependence case), $\varpi_n = n^{q/2-1- \varsigma q}$. The moments of the envelopes are assumed to be finite in \hyperref[C6]{(C6)}.}

{

}

The results in Theorem (ref) imply the asymptotic normality of the proposed estimator by Algorithm \hyperref[algo1]{1} and \hyperref[algo2]{2} by applying central limit theorems and Gaussian Approximation.

Consider the vector $\widetilde\zeta_t \stackrel{\mathrm{def}}{=} \operatorname{vec}\{(\zeta_{jk,t})_{(j,k)\in G}\}$, $\zeta_{jk,t} \stackrel{\mathrm{def}}{=} -\sigma^{-1}_{jk}\phi_{j,k}^{-1}\psi_{jk,t}^0$, and define the aggregated dependence adjusted norm as follows:

where $q\geqslant1$, and $\varsigma>0$. Moreover, define the following quantities

Define $L_1^\zeta = \{\Phi^\zeta_{2,\varsigma}\Phi^\zeta_{2,0} (\log|G|)^2\}^{1/\varsigma}$, $W_1^\zeta = \{(\Phi^\zeta_{3,0})^6+ (\Phi^\zeta_{4,0})^4\}\{\log(|G|n)\}^7$, $W_2^\zeta = (\Phi^\zeta_{2,\varsigma})^2\{\log(|G|n)\}^4$, $W_3^\zeta = [n^{-\varsigma} \{\log (|G|n)\}^{3/2} \Theta^\zeta_{q,\varsigma}]^{1/(1/2-\varsigma-1/q)}$, $N_1^\zeta =(n/\log|G|)^{q/2} (\Theta^\zeta_{q, \varsigma})^{q}$, $N_2^\zeta=n(\log|G|)^{-2}(\Phi^\zeta_{2,\varsigma})^{-2}$, $N_3^\zeta = \{n^{1/2}(\log|G|)^{-1/2}(\Theta^{\zeta}_{q, \varsigma}\})^{1/(1/2-\varsigma)}$.

Following Theorem (ref), a joint confidence region and the corresponding confidence interval for each component can be constructed via a block bootstrap method. In particular, the bootstrap statistics are defined by $\frac{1}{\sqrt{n}} \sum_{i=1}^{l_n} e_{j,i}\sum_{l=(i-1)b_n+1}^{ib_n} \widehat\zeta_{jk,l}$, where $e_{j,i}$'s are independent and identically distributed draws of standard normal random variables and are independent with respect to the data sample $(Z_{j,t})_{j=1}^J$. Recall that $\widehat\zeta_{jk,t}$ are pre-estimators with a certain range of accuracy. {More details can be found in Comment (ref) in the supplementary material.}

{

}

Simulation Study

In this section, we illustrate the performance of our proposed methodology under different simulation scenarios. The first part concerns the performance of the jointly selected penalty level over equations, and the second part discusses the simultaneous inference.

Estimation with a Jointly Selected Penalty Level

Consider the system of regression equations:

where $X_{t}\in{\rm I\!R}^K$. We generate $X_{t}$ independently from $\operatorname{N}(0,\Sigma)$, where $\Sigma_{k_1,k_2}=\gamma^{|k_1-k_2|}$, $\gamma=0.5$, $\varepsilon_{j,t}\stackrel{\operatorname{i.i.d.}}{\sim}\operatorname{N}(0,1)$. The coefficient vectors $\beta_j$ are assumed to be sparse. In particular, we divide the indices $\{1,\ldots,K\}$ evenly into blocks with fixed block size 5. $\beta_{jk}^0=10$ if $k$ and $j$ belong to the same block and 0 otherwise.

We take $n=100$, $\#$ of bootstrap replications = 5000. We set $J,K=50, 100$ and $150$. The prediction norm $| \widehat \beta_j - \beta^0_{j}|_{j,pr}$ and the Euclidean norm $| \widehat \beta_j - \beta^0_{j}|_{2}$ ratios are presented in Table (ref). The ratios measure the relative difference between the results using the penalty level determined from the equation-by-equation case and from the joint equation case ($\lambda_j$ and $\lambda$ are selected by the multiplier bootstrap procedure). In particular, a ratio smaller than $1$ indicates a better performance of using the jointly selected penalty level.

It is evident from Table (ref) that the proposed estimation procedure delivers much better performance in terms of the two measures. In particular, the superiority tends to be more evident (more than $10 \% $) with higher dimension of the covariates and more equations.

Still consider the system of regression equations as in (ref), but here we generate the data with dependency by following the Appendix D in ZW15_sup. In particular, assume the linear process such that $X_t=\sum_{\ell=0}^{\infty}A_\ell\xi_{t-\ell}$, with $A_\ell=(\ell+1)^{-\rho-1}M_\ell$, where $M_\ell$ are independently drawn from Ginibre matrices, i.e. all the entries of $M_\ell$ are i.i.d. $\operatorname{N}(0,1)$, and in practice the sum is truncated to $\sum_{\ell=0}^{1000}$. We set $\rho$ to be 1.0 for the weaker dependence and 0.1 for the stronger dependence cases respectively. Let $\xi_{k,t}=e_{k,t}(0.8e_{k,t-1}^2+0.2)^{1/2}$ where $e_{k,t}$ are i.i.d. distributed as $t(d)/\sqrt{d/(d-2)}$ and $t(d)$ is the Student's $t$ with degree of freedom $d$ (take $d=8$ for example). $\varepsilon_t$ are generated by following the same fashion independently.

We take $n=100$, $\#$ of bootstrap replications = 5000, $J,K=50, 100$ and $150$. {Based on bias-variance trade-off, several approaches were suggested to determine the optimal choice of $b_n$ for univariate case. Concerning the high-dimensional case, we propose to take the one which gives the lowest prediction norm as the optimal choice. Below we report the average prediction norm $J^{-1}\sum_{j=1}^J|\widehat \beta_j - \beta^0_{j}|_{j,pr}$ with several block sizes $b_n$ under different settings and the minimal ones are in bold. }

{From Table (ref), it is apparent that a larger block size is required for the stronger dependency case. Moreover, the choice also depends on the dimensionality, which is more evident for relatively weaker dependent data. We note that when $J=K=50$, $\rho=1.0$ the ordinary multiplier bootstrap (with $b_n=1$) produces 2.1003 as the average prediction norm, therefore we suggest $b_n=2$ for this case.}

The prediction norm $| \widehat \beta_j - \beta^0_{j}|_{j,pr}$ and the Euclidean norm $| \widehat \beta_j - \beta^0_{j}|_{2}$ ratios {(using the optimal $b_n$ suggested in Table (ref) for each case correspondingly)} are presented in Table (ref). Again we report the results with the jointly estimated $\lambda$ (selected by the algorithm proposed in section 3.2 based on multiplier block bootstrap) relative to using the single equation $\lambda_j$'s.

The results show that the coefficient estimation performance measured by both the prediction norm and the Euclidean norm is in favor of the joint penalty level with multiplier block bootstrap approach. The results are robust over different dimension cases with stronger or weaker dependency.

Simultaneous Inference

In this subsection we consider the following regression model for the purpose of simultaneous inference on the parameters within a system of equations

where $\alpha^0_j=\alpha^0$ for all $j$. Also, $\beta^0_j,\theta^0_j\in{\rm I\!R}^K$ are assumed to be sparse. In particular, we divide the indices $1,\ldots,K$ evenly into blocks with a fixed block size $5$, $\beta_{jk}^0$ and $\theta_{jk}^0$ are independently drawn from $\operatorname{Unif}[0,5]$ and $\operatorname{Unif}[0,0.25]$ respectively, if $k$ and $j$ belong to the same block and $0$ otherwise. The way to generate $X_t$, $\varepsilon_{t}$ and $v_t$ is same as the dependent data setting above.

We consider the sample size $n=100$. Our goal is to estimate and make inferences on the target variables $d_{j,t}$'s based on the procedure proposed in Section (ref). We evaluate and compare the empirical power and size performance of the confidence intervals constructed by the asymptotic distribution theory (ref), block bootstrap (ref) and the simultaneous confidence regions via block bootstrap (ref). The bootstrap statistics are computed based on 5000 replications and we also take the optimal block size according to the numerical comparison conducted above. Note that the case of $\alpha^0=0$ gives the size performance under the null hypothesis, while $\alpha^0$ uniformly lies in $[0,2.5]$ and $[0,5]$ illustrate the power results.

Table (ref) shows the average rejection rate of $H_0^j: \alpha_j^0=0$ over $j$ for individual (or multiple) inference and the rejection rate of $H_0: \alpha_1^0=\cdots=\alpha_J^0=0$ for simultaneous inference under different settings of $J,K$ and $\rho$. {Multiple testing procedure via step-down method (see e.g. romano2005exact,CCK13AoS), is considered to control the false positives in evaluating the power performance.} The rejection rates are computed over $1000$ simulation samples.

{It is shown that for individual inference our proposed individual bootstrap approach provides a closer size control to the nominal $\alpha$ and more powerful empirical rejection probabilities compared to constructing the confidence intervals by asymptotic normality in most of the cases. Moreover, the simultaneous inference outperforms the individual inference in size accuracy and in terms of the power performance, the multiple testing is relatively conservative after controlling the false positives. Overall, we observe that the results using bootstrap approach are robust over different dimension settings under either stronger or weaker dependency cases. }

Empirical Analysis: Textual Sentiment Spillover Effects

Financial markets are driven by information, and this is a well-known phenomenon among investors. More frequent news and availability of sentiment data allows study of the impact of firm-specific investor sentiment on market behavior such as stock returns, volatility and liquidity; see baker2006, tetlock2007, among others. Moreover, powerful statistical tools (e.g. LASSO-type estimators) are being used to model complex relationships among individuals. For example, audrino_sentiment analyze the influence of news on US and European companies by constructing a sparse predictive network via adaptive LASSO and related testing procedures. In this section the developed technology is applied to study textual sentiment spillover effects across individual stocks. This is different from the "equation-by-equation" analysis in audrino_sentiment, since we build up a system of regression equations and implement the estimation and the inference of the network jointly.

Data Source

The empirical study in this paper is carried out based on the financial news articles published on the NASDAQ community platform from January 2, 2015 to December 29, 2015 (252 trading days). The data were gathered via a self-written web scraper to automate the downloading process. The dataset is available at the Research Data Centre (RDC), Humboldt-Universit\"at zu Berlin. Moreover, unsupervised learning approaches are employed to extract sentiment variables from the articles. Two sentiment dictionaries: the BL option lexicon BL2004text and the LM financial sentiment dictionary LM2011text were used in zhang2016text. For each article $i$ (published on day $t$), the average proportion of positive/negative words using BL or LM lexica - $Pos_{j,i,t}^\text{BL}$, $Neg_{j,i,t}^\text{BL}$, $Pos_{j,i,t}^\text{LM}$, $Neg_{j,i,t}^\text{LM}$ - are considered as the text sentiment variables. Furthermore, the bullishness indicator for stock $j$ on day $t$ with the related articles $i=1,\ldots,m$ (based on a particular lexicon) is constructed by following af2004text

We refer to zhang2016text for more details about the data gathering and processing procedure. 63 individual stocks which are S&P 500 component stocks from 9 Global Industrial Classification Standard (GICS) sectors are considered. They are traded at NSDAQ Stock Exchange or NYSE. The list of the stock symbols and the corresponding company names can be found in Table (ref) in Appendix (ref) in the supplementary materials.

The daily log returns $R_{j,t}$ and log volatilities $\log(\sigma_{j,t}^2)$ for the stocks over the same time span are taken as response variables. More precisely, the gk1980vola range-based measure to represent the volatility level is employed:

where $u_{j,t} = \log(P_{j,t}^H) - \log(P_{j,t}^O), d_{j,t} = \log(P_{j,t}^L) - \log(P_{j,t}^O), r_{j,t} = \log(P_{j,t}^C) - \log(P_{j,t}^O)$, with $P_{j,t}^H, P_{j,t}^L, $, $P_{j,t}^O$, and $P_{j,t}^C$ denote the highest, lowest, opening and closing prices, respectively. In addition, the S&P 500 index returns and Chicago Board Options Exchange volatility index (VIX) are included as the state variables. The financial time series data were originally obtained from Datastream, and GICS sector information was found at Compustat.

Model Setting and Results

We now construct a network model to detect the spillover effects from sentiment variables to financial variables by

where $j=1,\ldots,J$ indicate the stock symbols, $B_{t}=(B_{1,t},\ldots,B_{J,t})^\top$ and $z_t$ includes the state variables.

It is of interest to make inferences on the parameters $\beta_j\in {\rm I\!R}^J$, $j=1,\ldots J$. Following the framework introduced in Section (ref), an estimation procedure with three steps needs to be implemented.

If for stock $j$, the sentiment variable of firm $k$ is selected into the active set after the individual significance test i.e., the null hypothesis $H_0^{jk}: \beta_{jk}=0$ is rejected under the block multiplier bootstrap procedure {(as discussed in Section (ref) we pre-determine $b_n=5$ by choosing the one gives the lowest prediction norm in the LASSO estimation in S1 on a grid search)}, then we put a directional edge from $k$ to $j$. As a result, we achieve a $0-1$ adjacency matrix describing the dependency network from sentiment variable to financial variable. Note that the diagonal elements in the matrix show the self-effect of stocks.

The graphical network for stock returns and volatility modelled by (ref) based on BL and LM lexica (from 01/02/15 to 12/29/15) is depicted in Figures (ref)-(ref).

Figures (ref)-(ref) depict the dependency networks among individual stocks. Given that the time series of returns and volatility are scaled and centered before implementing the estimation procedure, we find even denser spillover effects in the volatility analysis. This indicates the stock volatility is more sensitive to sentiment than returns. Moreover, the relationships between sectors are also of interest. The simultaneous confidence region constructed via the bootstrap approach introduced in Section (ref) may help us to detect whether the sentiment information from one sector has joint influence on the returns of the stocks in another sector. In particular, we look at the null hypothesis: $H_0^{S_1,S_2}:\, \beta_{jk}=0,\,\,\forall j\in S_1,\,k\in S_2$, where $S_1$ and $S_2$ represent two groups of stocks that belong to two sectors, respectively. The conclusion that the sentiment from sector $S_2$ has a joint effect on the returns or volatility of sector $S_1$ can be drawn if the null hypothesis is rejected with the simultaneous confidence region (ref) under the significance level = 0.05.

Figure (ref) describes the spillover effect network from sentiment to financial variables on the sector levels. In particular, the connections from energy to health care is found to be significant in the analysis of stock returns; while if volatility is focused on then the spillover effects from financials to health care, from information technology to energy, also from consumer discretionary to utilities are detected.

{

}

\vskip 2em \centerline{ \bf Supplementary Material} \vskip -1em \setcounter{subsection}{0} \vskip 2em