Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
127,626 characters · 15 sections · 84 citation commands
Empirical likelihood approach for high-dimensional moment restrictions with dependent data
\bibpunct{(}{)}{,}{a}{;}
\if11 {
\affil[1]{\it Joint Laboratory of Data Science and Business Intelligence, Southwestern University of Finance and Economics, Chengdu, Sichuan, China} \affil[2]{\it Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, China} \affil[3]{\it Department of Economics, Chinese University of Hong Kong, Hong Kong SAR, China}
} \fi
\if01 {
} \fi
\onehalfspacing
{\sl keywords}: $\alpha$-mixing, asymptotic analysis, confidence region, high dimensionality, penalized likelihood
Econometrics began almost a hundred years ago with a focus on deciphering business cycles and economic fluctuations, as demonstrated by early researchers fisher1925our, frisch1933propagation, and tinbergen1939. These pioneers carried out their quantitative analyses with limited macro-level, low-frequency datasets. However, with the progress in information technology, economists now have access to extensive long-term data series that reflect the state of real economies and financial markets from diverse angles. Despite this, the complexity of the economic landscape continues to escalate, characterized over the last century by exceptional growth interspersed with recessions, wars, and crises. The age-old inquiry into the interplay between economic variables and their evolving paths still holds a prominent place in the current big data era.
An important class of economic models is characterized by moment conditions. To fix ideas, let $\{{\mathbf x}_t\}_{t=1}^n$ be $d$-dimensional random vectors, and $\boldsymbol{\theta}=(\theta_1,\ldots,\theta_p)^{{ \mathrm{\scriptscriptstyle \top} }}$ be a $p$-dimensional parameter in a parameter space $\boldsymbol{\Theta}$. For an $r$-dimensional estimating function ${\mathbf g}({\mathbf x};\boldsymbol{\theta})=\{g_1({\mathbf x};\boldsymbol{\theta}),\dots,g_r({\mathbf x};\boldsymbol{\theta})\}^{ \mathrm{\scriptscriptstyle \top} }$, the information for the true parameter $\boldsymbol{\theta}_0$ is identified by the moment condition
for each $t=1,\ldots,n$. In the traditional asymptotic analysis, it is assumed that the sample size $n$ approaches infinity while the model remains unchanged. Economic theory, however, often relies on the orthogonality of variables as a justification for moment conditions, but it generally lacks clarity regarding which variables should be included or which moment conditions should be utilized. For example, angrist1991does and eaton2011 involve a large number of moments, and blundell1993we and fan2014endogeneity further contain many parameters. Beyond these applications in labor economics, international trade, and household consumption, multivariate time series models provide fertile ground for the development of complex economic models.
This study is driven by three commonly utilized multivariate time series models, as detailed in lutkepohl2005new's monograph. First, the vector autoregressive (VAR) model sims1980macroeconomics is a multivariate extension of the univariate autoregressive (AR) model. In a VAR model composed of $d$ variables, the total number of slope coefficients is $p= O(d^2)$. This is referenced as our Example 1 in Section (ref). Conventionally, the VAR model only includes a limited number of variables, which is not designed to handle a modern large system of variables.
The second example is concerning the impulse response function (IRF), which describes the effect of an external shock to a target variable. IRF is a key object of macroeconomic interest; see the survey by nakamura2018identification. Traditionally, IRF is implied by a fully and “correctly specified” VAR model. Recently, Jorda2005 introduced local projection (LP), a much simpler method which involves a series of single-equation regressions over the forecast horizons $h = 0, 1, \ldots, H$. If each individual regression includes $d$ regressors, the system yields $p = d(H+1)$ slope coefficients to be estimated.
In addition to the above two examples of conditional mean models, volatility is fundamental for financial risk management. The most recognized univariate time series volatility models are the autoregressive conditional heteroskedasticity (ARCH) model engle1982autoregressive and its generalized version (GARCH) bollerslev1986generalized. In typical financial markets, numerous assets are traded daily. To handle this complexity, multivariate volatility models, such as the multivariate ARCH (MARCH) and multivariate GARCH (MGARCH) EngleandKroner1995, have been developed. As will be explained in Example 3 in Section (ref), these models require high-dimensional coefficients to capture the interdependencies among various assets.
Moment restrictions serve as a comprehensive framework for identifying the parameters of interest across all the above-mentioned three multivariate time series models. Within the context of these moment constraints, the generalized method of moments (GMM) is widely regarded as the default econometric estimation approach hansen1982large, hansen1982generalized. Although Hansen's two-step GMM achieves asymptotic efficiency under the traditional asymptotic scenario with fixed $r$ and $p$, its performance in finite samples can be unsatisfactory, as is highlighted by altonji1996small. This issue primarily arises from the inversion of the estimated covariance matrix of the estimating functions. Recently, cheng2023weight investigate the bias of GMM in cases where $p$ is proportional to, but smaller than $n$. When $p$ exceeds $n$, belloni2018high propose using a sup-norm objective criterion function to manage many moments, departing from GMM’s quadratic form.
An alternative to GMM is empirical likelihood (EL) Owen(1988), QinandLawless(1994), which can be viewed through the lens of information theory kitamura1997information. Unlike methods that require explicit computation of the covariance matrix of the estimating functions, EL benefits from reduced variance due to higher-order enhancements newey2003higher. When dealing with models incorporating many moments, the EL criterion function can be augmented with a penalized approach to regularize both the multitude of moments and high-dimensional coefficients. This branch of theoretical properties, specifically for independently and identically distributed (i.i.d.) observations, has been extended by otsu2007penalized, leng2012penalized, shi2016econometric, and CTW_2018. Penalized method theory is relatively straightforward with i.i.d. data. To our knowledge, this paper is the first to explore the EL methodology in the context of high-dimensional temporally dependent data.
In classical low-dimensional settings, kitamura1997empirical advocates using blocks of time series observations in EL to preserve the temporal dependence and recover the Wilks' phenomenon. This blocking technique is employed by CCC(2015) under a moderately high-dimensional environment where $r/n\to 0$. However, the blocking technique is inconvenient when dealing with high-dimensionality in both the parameters and moments. This paper maintains the simplest approach: treating the EL as if the data is i.i.d., which can be interpreted as the marginal EL for time series. Though marginal likelihood estimation has been used in low-dimensional models levenbach1972estimation and has been applied in financial applications (stambaugh1997analyzing, patton2006estimation) for integration of varying lengths, this paper is the first to investigate the scheme of marginalization in the framework of EL for time series, distinct from the approaches by ChangTangWu2013,chang2016local.
We establish rigorous asymptotic theory for the penalized EL (PEL) with many parameters and many moments under time dependence of $\alpha$-mixing. To accommodate stronger temporal dependence, we introduce (in Condition (ref)) a diverging quantity $L_n$ within the $\alpha$-mixing coefficient. To address the high-dimensional challenges in this framework, we bound the tail probabilities of certain key statistics by novel inequalities, which are constructed via self-normalized sum inequalities JingShaoWang2003. We show that under sparsity conditions, the PEL approach delivers consistent estimation, and the resulting PEL estimator is asymptotically normally distributed. For investigations focused on a low-dimensional parameter, we further employ a projected PEL (PPEL) method Chang2020 to eliminate the bias induced by the high-dimensional nuisance parameters, thereby reinstating the usual inference method using the $t$-statistic.
This procedure enables the application of the method to a wide range of multivariate time series models with many parameters, including VAR, LP, and MARCH/MGARCH models as discussed. Extensive Monte Carlo simulations show that our method performs well in finite sample. We employ it to further study the persistence and spillover of the USA's sectoral inflation, the magnitude of the fiscal multiplier, and the volatility network of China's banking industry.
This paper stands on the large literature of time series, empirical likelihood, and high-dimensional estimation. It leverages the penalized estimation in CTW_2018 and Chang2020, developed under the i.i.d. setting. These steps are carried over and adapted to the environment with temporal dependence. Compared to the GMM-type alternative, belloni2018high is developed in the i.i.d. environment, which is critical for their sample splitting; sample splitting in time series is much more challenging and may adversely affect finite sample performance when the time length is moderate in practice. On the other hand, the literature of high-dimensional time series has witnessed specific proposals for standalone models. For example, shi2022 and mei2024lasso focus on high-dimensional time series dense regressions and sparse regressions, respectively. Under the assumption of Gaussian errors, kock2015oracle develop Lasso-type penalized estimation for VAR models. caner2018high attack high-dimensional GMM-type models, where their linear setting facilitates the estimation of the large weighting matrix. adamek2024local provide theory for a single-equation regression with many covariates for local projection, and deal with ridge-type regularization in fixed dimension. Our procedure provides a unified framework to handle models defined by moment conditions.
The rest of the paper is organized as follows. Section (ref) sets up the model and the technical conditions. The consistency and asymptotic distribution of the PEL estimator are established in Section (ref), and the asymptotic normality of PPEL is presented in Section (ref). We carry out Monte Carlo simulations in Section (ref), and showcase our method in three empirical applications in Section (ref). The code and data are available at the GitHub repository: \url{https://github.com/JinyuanChang-Lab/PenalizedELwithDependentData}. All proofs are relegated to the Appendices.
{\bf Notation}. We use the abbreviations “w.p.a.1” and “w.r.t” to denote, respectively, with probability approaching one and with respect to. For any real number $x$, define $ \lfloor x \rfloor = \max \{ q \in \mathbb{Z} : q \leq x \} $, where $\mathbb{Z}$ denotes the set of all integers. For two sequences of positive numbers $\{a_n\}$ and $\{b_n\}$, we write $a_n\lesssim b_n$ or $b_n\gtrsim a_n$ if there exists a positive constant $c$ such that $\lim\sup_{n\to \infty} a_n/b_n\leq c$, and $a_n\asymp b_n$ if and only if $a_n\lesssim b_n$ and $b_n\lesssim a_n$ hold simultaneously. We write $a_n\ll b_n$ or $b_n\gg a_n$ if $\lim\sup_{n\to \infty} a_n/b_n=0$. Let “vec” and “vech” be the vector operators that stack the columns of a matrix and the upper triangular part of a matrix, respectively, into a vector. For a positive integer $q$, we write $[q]=\{1,\ldots,q\}$, and let ${\mathbf I}_{q}$ be the $q\times q$ identity matrix. For a $q\times q$ symmetric matrix ${\mathbf Q}$, denote by $\lambda_{\min}({\mathbf Q})$ and $\lambda_{\max}({\mathbf Q})$ the smallest and largest eigenvalues of ${\mathbf Q}$, respectively. For a $q_1\times q_2$ matrix ${\mathbf B}=(b_{i,j})_{q_1\times q_2}$, let ${\mathbf B}^{{ \mathrm{\scriptscriptstyle \top} }}$ be its transpose, $|{\mathbf B}|_{\infty}=\max_{i\in[q_1],j\in[q_2]}|b_{i,j}|$ be the sup-norm, and $\|{\mathbf B}\|_2=\lambda_{\max}^{1/2}({\mathbf B}^{\otimes2})$ be the spectral norm with ${\mathbf B}^{\otimes2}={\mathbf B}{\mathbf B}^{{ \mathrm{\scriptscriptstyle \top} }}$. Specifically, when $q=1$, we use $|{\mathbf B}|_{1}=\sum_{i=1}^{q_1}|b_{i,1}|$ and $|{\mathbf B}|_2=(\sum_{i=1}^{q_1}b_{i,1}^2)^{1/2}$ to denote the $L_1$-norm and $L_2$-norm of the vector ${\mathbf B}$. For two square matrices ${\mathbf Q}_1$ and ${\mathbf Q}_2$, we say ${\mathbf Q}_1\leq {\mathbf Q}_2$ if $({\mathbf Q}_2-{\mathbf Q}_1)$ is a positive semi-definite matrix. The population mean is denoted by $\mathbb{E}(\cdot)$, and the sample mean is denoted by $\mathbb{E}_n(\cdot)=n^{-1}\sum_{t=1}^{n}(\cdot)$. For a given index set $\mathcal{L}$, let $|\mathcal{L}|$ be its cardinality. For a generic multivariate function ${\mathbf h}(\cdot;\cdot)$, we denote by ${\mathbf h}_{{ \mathcal{\scriptscriptstyle L} }}(\cdot;\cdot)$ the subvector of ${\mathbf h}(\cdot;\cdot)$ collecting the components indexed by $\mathcal{L}$. Analogously, we write ${\mathbf a}_{\mathcal{L}}$ as the corresponding subvector of a vector ${\mathbf a}$. For simplicity and when no confusion arises, we use the generic notation ${\mathbf h}_t(\boldsymbol{\theta})$ as equivalent to ${\mathbf h}({\mathbf x}_t;\boldsymbol{\theta})$, and $\nabla_{\boldsymbol{\theta}}{\mathbf h}_t(\boldsymbol{\theta})$ for the first-order partial derivative of ${\mathbf h}_t(\boldsymbol{\theta})$ w.r.t $\boldsymbol{\theta}$. Denote by $h_{t,k}(\boldsymbol{\theta})$ the $k$-th component of ${\mathbf h}_t(\boldsymbol{\theta})$. Let $\bar{{\mathbf h}}(\boldsymbol{\theta})=\mathbb{E}_n\{{\mathbf h}_t(\boldsymbol{\theta})\}$, and write its $k$-th component as $\bar{h}_k(\boldsymbol{\theta})=\mathbb{E}_n\{h_{t,k}(\boldsymbol{\theta})\}$. Analogously, let ${\mathbf h}_{t,{ \mathcal{\scriptscriptstyle L} }}(\boldsymbol{\theta})={\mathbf h}_{{ \mathcal{\scriptscriptstyle L} }}({\mathbf x}_t;\boldsymbol{\theta})$ and $\bar{{\mathbf h}}_{{ \mathcal{\scriptscriptstyle L} }}(\boldsymbol{\theta})=\mathbb{E}_n\{{\mathbf h}_{t,{ \mathcal{\scriptscriptstyle L} }}(\boldsymbol{\theta})\}$.
In this paper, we build up our theory with $\alpha$-mixing time dependence. Let $\mathcal{F}_{-\infty}^u$ and $\mathcal{F}_{u}^{\infty}$ be the $\sigma$-fields generated by $\{{\mathbf x}_{t}\}_{t\leq u}$ and $\{{\mathbf x}_{t}\}_{t\geq u}$, respectively. The $\alpha$-mixing coefficient of the sequence $\{{\mathbf x}_t\}$ at lag $k$ is defined as
for each $k\geq 1$. The notion of $\alpha$-mixing in $\eqref{eq:alpha-mixing}$ broadly characterizes serial dependence. Specifically, we impose the following assumption as in chang2024optimal in our study.
Condition (ref) does not require $\{{\mathbf x}_t\}_{t=1}^{n}$ to be strictly stationary. For an independent sequence $\{{\mathbf x}_t\}_{t=1}^{n}$, we can select $L_n=1/2$ and $\varphi=\infty$. For an $L_n$-dependent sequence $\{{\mathbf x}_t\}_{t=1}^{n}$, we can select $\varphi=\infty$. A variety of time series models that are routinely used in economics and finance are covered by Condition (ref). For example, under some regularity conditions, the autoregressive-moving-average (ARMA) processes, the stationary Markov chains FanYao_2003, and the stationary GARCH models CarrascoChen_2002 satisfy $\alpha$-mixing with the exponentially decaying coefficient ($L_n=1$ and $\varphi=1$). This condition further covers their multivariate generalizations of VAR and MGARCH (MGARCH includes MARCH as a special case); see HP09, BFS11 and Wong2020.
The quantity $L_n$ involved in Condition (ref) accommodates practical scenarios of big data collected over time. Write ${\mathbf x}_t=(x_{t,1},\ldots,x_{t,d})^{{ \mathrm{\scriptscriptstyle \top} }}$. Consider the simple case when each univariate time series $\{x_{t,i}\}_{t=1}^n$ for $i \in [d]$ is $\alpha$-mixing with exponentially decaying $\alpha$-mixing coefficients, while those $d$ sequences are mutually independent. Theorem 5.1 of bradley2005basic indicates that $\alpha_n(k)$ defined in (ref) satisfies $\alpha_n(k)\leq d\exp(-ck)$ for some universal constant $c>0$, which implies Condition (ref) holds for $\varphi=1$ and $L_n\asymp \log d$. Our novel framework also covers high-frequency time series models. Suppose the observed vector process is generated from ${\mathbf x}_{t}={\mathbf P}{\mathbf z}_{t\delta}$, for $t\in[n]$, where ${\mathbf P}\in\mathbb{R}^{ d\times q}$ is a loading matrix and $\delta>0$ is the sampling interval. Let the latent vector process ${\mathbf z}_s=(z_{s,1},\ldots,z_{s,q})^{{ \mathrm{\scriptscriptstyle \top} }}$ consist of $q$ independent processes, where each $z_{s,i}$ for $i\in[q]$ follows the diffusion model ${\rm d} z_{s,i} = \tilde{\mu}_i(z_{s,i})\,{\rm d}s+\tilde{\sigma}_i(z_{s,i})\,{\rm d}{W}_{s,i}$, with a univariate standard Brownian motion ${W}_{s,i}$ and two parametric functions $\tilde{\mu}_i(\cdot)$ and $\tilde{\sigma}_i(\cdot)$. When $\tilde{\mu}_i(\cdot)$ and $\tilde{\sigma}_i(\cdot)$ satisfy certain conditions as those in Lemma 4 of ait2004estimators, the observed $\{{\mathbf x}_t\}_{t=1}^{n}$ satisfies Condition (ref) with $\varphi=1$ and $L_n=\delta^{-1}$, where $L_n$ diverges if $\delta\rightarrow 0$ as $n\rightarrow\infty$.
We first set up the estimation procedure. We are interested in the $p$-dimensional parameter $\boldsymbol{\theta}_{0}\in \boldsymbol{\Theta}$ defined as the solution of the $r$ moment conditions (ref). Based on Owen(1988), Owen(1990)'s seminal idea, given the estimating equations $\{{\mathbf g}_t(\cdot)\}_{t=1}^n$, QinandLawless(1994) define the EL as
Maximizing $L(\boldsymbol{\theta})$ can be equivalently carried out via the corresponding dual problem, and its optimizer is the EL estimator:
where $\boldsymbol{\lambda}=(\lambda_1,\ldots,\lambda_r)^{{ \mathrm{\scriptscriptstyle \top} }}$ and $\hat{\Lambda}_n(\boldsymbol{\theta})=\{\boldsymbol{\lambda}\in\mathbb{R}^r:\boldsymbol{\lambda}^{{ \mathrm{\scriptscriptstyle \top} }}{\mathbf g}_t(\boldsymbol{\theta})\in \mathcal{V} \text{ for any }t\in[n]\}$ for an open interval $\mathcal{V}$ containing zero.
Most economic and financial time series typically consist of a few hundred observations or more. Meanwhile, models that characterize economic interactions are inherently complex and high-dimensional. Regularization is crucial for accurately estimating many parameters with limited time datasets. Shrinkage serves as a useful statistical technique for dimension reduction, with its effectiveness depending on the nature of the data and models. Arguably, sparsity is an extensively used assumption in high-dimensional models, exemplified by the success of Lasso tibshirani1996regression, SCAD FanLi2001 and MCP Zhang2010, which have been applied across various scientific fields.
Consider a model with the number of estimating equations $r$ and the number of parameters $p$ both potentially larger than the sample size $n$. Such a high-dimensional model can be estimated by CTW_2018's PEL method:
where two penalty functions $P_{1,\pi}(\cdot)$ and $P_{2,\nu}(\cdot)$ with tuning parameters $\pi$ and $\nu$ are appended to the dual problem (ref). For any penalty function $P_\tau(\cdot)$ with a tuning parameter $\tau$, let $\rho(t;\tau)=\tau^{-1}P_{\tau}(t)$ for any $t \in [0,\infty)$ and $\tau \in (0,\infty)$. Assume the penalty functions $P_{1,\pi}(\cdot)$ and $P_{2,\nu}(\cdot)$ belong to the following class as in LvandFan(2009):
The PEL estimator (ref) is formulated as if ${\mathbf x}_t$ is i.i.d., which contrasts with kitamura1997empirical: “... studies the method of empirical likelihood in models with weakly dependent processes. In such cases, if the likelihood function is formulated as if the data process were independent, obviously empirical likelihood fails.” To restore the Wilks' phenomenon, kitamura1997empirical proposes the blocking technique for low-dimensional EL estimation under a fixed $r$, and this method is adopted by CCC(2015) under $r\to \infty $ with $r/n\to 0$. However, blocking with a length $b$ reduces the effective sample size from $n$ to $\lfloor n/b \rfloor$. Asymptotic theory for a long vector of ${\mathbf x}_t$ would request a large block size to cope with the variable in ${\mathbf x}_t$ of the maximum temporal dependence. In the finite sample, a big $b$ for all time series would substantially reduce the effective sample size; on the other hand, choosing a block size for each time series would involve many more additional tuning parameters. Block preserves nice statistical properties in low-dimensional time series models, but it is inconvenient in high-dimensional contexts. Our marginal EL in (ref) circumvents the choice of block sizes. We will maintain the PEL formulation as for the i.i.d. data and develop the theory accordingly.
When economic theory provides no clear guidance about which variables are the most relevant or which moments are the most informative, shrinkage methods are helpful as a data-driven device for variable and moment selection. Given the $p$-dimensional true parameter $\boldsymbol{\theta}_0 =(\theta_{0,1},\dots,\theta_{0,p})^{{ \mathrm{\scriptscriptstyle \top} }} $, let $\mathcal{S} = \{k\in[p]:\theta_{0,k}\neq 0 \}$ be the active set of cardinality $s=|\mathcal{S}|$, where $\mathcal{S}$ marks the location of the non-zero parameters with $s\ll p$. Before any attempt at estimation, the parameter of interest must be identifiable from the population model.
This assumption means that the expected values of the estimating functions at the true parameter $\boldsymbol{\theta}_0$ are significantly different from those outside a narrow vicinity of the active coefficient $\boldsymbol{\theta}_{0,{ \mathcal{\scriptscriptstyle S} }}$. The sup-norm on moment conditions is a basic necessity that allows for the inclusion of many weak or entirely irrelevant moments. Under this condition, the parameter is identified locally as described by chen2014local.
We move on to the proceeding conditions which are standard regularity assumptions in the literature. For any index set $\mathcal{F}\subset[r]$ and $\boldsymbol{\theta}\in\boldsymbol{\Theta}$, define $\widehat{{\mathbf V}}_{{ \mathcal{\scriptscriptstyle F} }}(\boldsymbol{\theta})=\mathbb{E}_n\{{\mathbf g}_{t,{ \mathcal{\scriptscriptstyle F} }}(\boldsymbol{\theta})^{\otimes2}\}$. When $\mathcal{F}=[r]$, we write $\widehat{{\mathbf V}}_{{ \mathcal{\scriptscriptstyle F} }}(\boldsymbol{\theta})=\widehat{{\mathbf V}}(\boldsymbol{\theta})$ for conciseness.
Condition (ref)(a) restricts to be finite the $(2\gamma)$-th population moments of the estimating functions, and the $\gamma$-th sample moments to be $O_{\rm p}(1)$ in a uniform manner. At the true parameter $\boldsymbol{\theta}_0$, Condition (ref)(b) requires non-degenerate long-run variance, and Condition (ref)(c) ensures a well-behaved covariance matrix $\mathbb{E}\{\widehat{{\mathbf V}}(\boldsymbol{\theta}_{0})\}$ for the estimating functions, which is satisfied if the eigenvalues of ${\rm Var}\{{\mathbf g}_1(\boldsymbol{\theta}_0)\},\ldots,{\rm Var}\{{\mathbf g}_n(\boldsymbol{\theta}_0)\}$ are uniformly bounded away from zero and infinity. Next, Condition (ref) mimics those in Condition (ref) by regularizing the derivatives of the estimating functions.
Finally, we deal with the choice of the penalty functions. Consider a generic penalty $P_{1,\pi}(\cdot)$ and let $a_n = \sum_{k=1}^{p}P_{1,\pi}(|\theta_{0,k}|)$. Define $b_n$ = $\max\{a_n,\nu^2\}$, which is the larger value between $a_n$ and the square of the tuning parameter $\nu$ attached to the second penalty function $P_{2,\nu}(\cdot)$. In order to control the shrinkage bias induced by $P_{1,\pi}(\cdot)$ on $\hat{\boldsymbol{\theta}}_n$, suppose there exist $\chi_n\to 0$ and $c_n\to 0$ with $c_n \gg b_n^{1/2} $ such that
Under the assumption $b_n\ll \min_{k\in\mathcal{S}}|\theta_{0,k}|^2$, we can replace (ref) by
for some constant $c\in(0,1)$. If we select $P_{1,\pi}(\cdot)$ as an asymptotically unbiased penalty such as SCAD or MCP, we have $\chi_n=0$ in (ref) when
To simplify the presentation, we assume that (ref) holds and $\chi_n=0$ in (ref). It provides the minimum signal level on the nonzero components in $\boldsymbol{\theta}_0$.
Regarding the second penalty term, we write $\rho_2 (t;\nu)=\nu^{-1}P_{2,\nu}(t)$ for $P_{2,\nu}(t) \in \mathcal{P}$. Since $\rho_2 '(0^{+};\nu)$ is independent of $\nu$, we denote $\rho_2 '(0^{+};\nu)$ by $\rho_2 '(0^{+})$ for simplicity. For any $\boldsymbol{\theta} \in \boldsymbol{\Theta}$, define
for some constant $C_* \in (0,1)$. Given $\nu$, the complexity of the moment conditions can be controlled by a “sparsity of moments” index $\ell_n$ such that
with $c_{n}\rightarrow 0$ satisfying $ c_{n} \gg b_n^{1/2} $, where $\aleph_n=n^{-3\varphi/(6\varphi+2)}(\log r)^{1/2}$.
The technical conditions in the previous section allow us to proceed with the asymptotic properties of PEL. To facilitate analysis under the listed assumptions, we refine the problem (ref) as
where $\boldsymbol{\Theta}_*=\{\boldsymbol{\theta}\in \boldsymbol{\Theta}: |\boldsymbol{\theta}_{{ \mathcal{\scriptscriptstyle S} }}-\boldsymbol{\theta}_{0,{ \mathcal{\scriptscriptstyle S} }}|_{\infty} \leq c_*,|\boldsymbol{\theta}_{{ \mathcal{\scriptscriptstyle S} }^{\rm c}}|_1\leq \aleph_n\}$ for a fixed $c_*>0$, where $\aleph_n=n^{-3\varphi/(6\varphi+2)}(\log r)^{1/2}$ has been defined in (ref). The refinement comes from $c_*$, which restricts the true value to be a point in the closed parameter space, and the sparsity is controlled by $\aleph_n$. Next, let $\hat{\boldsymbol{\lambda}} = (\hat{\lambda}_{1},\ldots,\hat{\lambda}_{r})^{ \mathrm{\scriptscriptstyle \top} }$ be the $r$-dimensional vector of Lagrange multiplier defined at $\hat{\boldsymbol{\theta}}_n${\rm:}
and $\mathcal{R}_{n}={\rm supp}(\hat{\boldsymbol{\lambda}})$ be its support. Define $\hat{\boldsymbol{\eta}}=(\hat{\eta}_1,\ldots,\hat{\eta}_{r})^{{ \mathrm{\scriptscriptstyle \top} }}$ with
It holds by the Karush-Kuhn-Tucker condition that $\eta_j=\nu\rho_2'(|\hat{\lambda}_j|;\nu)\mbox{\rm sgn}(\hat{\lambda}_j)$ for $\hat{\lambda}_{j}\neq 0$, and the subdifferential at $\hat{\lambda}_j$ must include the zero element Bertsekas(1997) for $\hat{\lambda}_j=0$. We impose the next condition.
Condition (ref)(a) is a technical assumption used to derive the convergence rate of the Lagrange multiplier $\hat{\boldsymbol{\lambda}}(\hat{\boldsymbol{\theta}}_n)$ associated with $\hat{\boldsymbol{\theta}}_n$; see the proof of Theorem (ref) in Appendix (ref). Condition (ref)(b) requires that w.p.a.1 the nonzero $\hat{\eta}_j$ does not lie on the boundary, which is satisfied by continuous random variables. Condition (ref) makes sure that $\hat{\boldsymbol{\lambda}}(\boldsymbol{\theta})$ is continuously differentiable at $\hat{\boldsymbol{\theta}}_n$ w.p.a.1; see Lemma (ref) in Appendix (ref).
The conditions up to this point are sufficient to guarantee the consistency of the PEL estimator.
Theorem (ref) provides the consistency of the PEL estimator: $\hat{\boldsymbol{\theta}}_{n}$ in the true active set $\mathcal{S}$ converges in probability to its true value at rate $b_n^{1/2}$, whereas the coefficients in the inactive set $\mathcal{S}^{\rm c}$ are shrunken to exactly zero w.p.a.1. The orders involved in the statement give the admissible range of the tuning parameters. Recall $\aleph_n=n^{-3\varphi/(6\varphi+2)}(\log r)^{1/2}$ and $b_n$ = $\max\{a_n,\nu^2\}$ with $a_n = \sum_{k=1}^{p}P_{1,\pi}(|\theta_{0,k}|)$. Due to $a_n \lesssim s\pi$, Theorem (ref) requires that the tuning parameters $(\nu,\pi)$ satisfy
with $\ell_n\ll\min\{s^{-1}n^{-2/\gamma} \aleph_n^{-1}, s^{-3/2} \aleph_n^{-1/2} \}$ and $s\ll n^{\delta}$ for some $\delta= \min\{3\varphi/(6\varphi+2)-2/\gamma, \varphi/(6\varphi+2)\}$. There is a tradeoff between $s$ and $r$. For example, if $s\asymp n^{\kappa}$ with $\kappa\in[0,\delta)$ is of polynomial order of $n$, Theorem (ref) ensures the consistency of $\hat{\boldsymbol{\theta}}_n$ even if the number of moments $r$ diverges exponentially fast under the rate $ \log r\ll \min \{ n^{3\varphi/(3\varphi+1)-4/\gamma-2\kappa}, n^{3\varphi/(3\varphi+1)-6\kappa}, n^{\varphi/(3\varphi+1)}L_n^{-\varphi} \} $ with $L_n\ll n^{1/(3\varphi+1)}(\log n)^{-1/\varphi}$.
While Theorem (ref) gives the rate of convergence, we further characterize the asymptotic distribution of $\hat{\boldsymbol{\theta}}_{n,{ \mathcal{\scriptscriptstyle S} }}$. For any index set $\mathcal{F}\subset[r]$ and $\boldsymbol{\theta}\in\boldsymbol{\Theta}$, denote the long-run covariance matrix of $\{{\mathbf g}_{t,{ \mathcal{\scriptscriptstyle F} }}(\boldsymbol{\theta})\}_{t=1}^{n}$ by $\boldsymbol{\Xi}_{{ \mathcal{\scriptscriptstyle F} }}(\boldsymbol{\theta})={\rm Var}\{n^{-1/2} \sum_{t=1}^{n}{\mathbf g}_{t,{ \mathcal{\scriptscriptstyle F} }}(\boldsymbol{\theta})\},$ and denote $\boldsymbol{\Xi}(\boldsymbol{\theta}):=\boldsymbol{\Xi}_{{ \mathcal{\scriptscriptstyle F} }}(\boldsymbol{\theta})$ when $\mathcal{F}=[r]$.
Condition (ref)(a) corresponds to the sparse Riesz condition (chen2008extended, zhang2008sparsity, chang2022culling) and it is adapted to our setting for handling high-dimensional dependent data. Condition (ref)(b) ensures that the eigenvalues of the long-run covariance matrix of the sequence $ \{{\mathbf g}_{t}(\boldsymbol{\theta}_0)\}_{t=1}^{n}$ are uniformly bounded away from zero and infinity.
Theorem (ref) below shows that the PEL estimator for the nonzero components of $\boldsymbol{\theta}_0$ is asymptotically normal, with an asymptotic bias that deviates from zero due to the high dimensionality. In addition to (ref), the admissible range of the tuning parameters $(\nu,\pi)$ in Theorem (ref) is slightly narrowed. Define
Here is the statement.
Though in principle the bias term $\hat{\boldsymbol{\psi}}_{{ \mathcal{\scriptscriptstyle R} }_n}$ is estimable, it involves approximation to multiple components that are difficult to formulate and compute. Theorem (ref) is a property of the PEL estimator for $\boldsymbol{\theta}_{0,{ \mathcal{\scriptscriptstyle S} }}$, and it is uninformative about the model parameters indexed in $\mathcal{S}^{\rm c}$. Consequently, we do not advocate using Theorem (ref) for statistical inference. Instead, if we are interested in testing low-dimensional components of the model parameter, which is commonly the use case in applied econometrics, we recommend PPEL in the following section.
We focus on the inference about $\boldsymbol{\theta}_{{ \mathcal{\scriptscriptstyle M} }}\in\mathbb{R}^m$ for some small subset $\mathcal{M}\subset[p]$ with $|\mathcal{M}|=m$. In economic and financial applications, most often the inference falls in a single parameter under which $m=1$. Hence, we assume $m$ is a fixed integer for simplicity; there is no technical difficulty in allowing it to diverge with the sample size $n$. We write $\boldsymbol{\theta}=(\boldsymbol{\theta}_{{ \mathcal{\scriptscriptstyle M} }}^{ \mathrm{\scriptscriptstyle \top} },\boldsymbol{\theta}_{{ \mathcal{\scriptscriptstyle M} }^{\rm c}}^{ \mathrm{\scriptscriptstyle \top} })^{ \mathrm{\scriptscriptstyle \top} }$, where $\boldsymbol{\theta}_{{ \mathcal{\scriptscriptstyle M} }}\in\mathbb{R}^m$ contains the low-dimensional components of interest, and $\boldsymbol{\theta}_{{ \mathcal{\scriptscriptstyle M} }^{\rm c}}\in\mathbb{R}^{p-m}$ contains the nuisance parameters. We consider an inferential procedure similar to that in chang2022culling. Let ${\mathbf A}_n=({\mathbf a}_1^n,\dots,{\mathbf a}_m^n)^{ \mathrm{\scriptscriptstyle \top} }\in\mathbb{R}^{m\times r}$, which is found row by row via the optimization problem
where $\varsigma$ is a tuning parameter, $\hat{\boldsymbol{\theta}}_n $ is the PEL estimator from ((ref)), and $\{\boldsymbol{\xi}_{k}\}_{k=1}^m$ consists of the canonical basis of the linear space $\mathcal{M}_{\boldsymbol{\xi}}=\{\mathbf{b}=(b_1,\dots,b_p)^{ \mathrm{\scriptscriptstyle \top} }:b_j=0 \text{ for any } j=m+1,\dots,p\}$, i.e., $\boldsymbol{\xi}_k$ is chosen such that its $k$-th component is $1$ and all other components are $0$. Define ${\mathbf f}^{{\mathbf A}_n}(\cdot;\cdot)={\mathbf A}_n{\mathbf g}(\cdot;\cdot)$, and then $ {\mathbf f}^{{\mathbf A}_n}(\cdot;\cdot)$ are the new $m$-dimensional estimating functions. By construction, the influence of the nuisance parameters is projected out by ${\mathbf a}^n_k$.
The low-dimensional moment functions ${\mathbf f}^{{\mathbf A}_n}(\cdot;\cdot)$ substantially reduce the dimension $r$ all the way to a much smaller $m$. The EL constructed with ${\mathbf f}^{{\mathbf A}_{n}}(\cdot;\cdot)$, instead of ${\mathbf g}(\cdot;\cdot)$, is used for the inference about $\boldsymbol{\theta}_{{ \mathcal{\scriptscriptstyle M} }}$. Specifically, let
where $\hat{\boldsymbol{\theta}}_{n,{{ \mathcal{\scriptscriptstyle M} }}^{\rm c}}$ is the nuisance component of $\hat{\boldsymbol{\theta}}_n$ from (ref). Maximizing ((ref)) yields the PPEL estimator $\tilde{\boldsymbol{\theta}}_{{ \mathcal{\scriptscriptstyle M} }}$, which can also be obtained by solving the corresponding dual problem
where $\widehat{\boldsymbol{\Theta}}_{{ \mathcal{\scriptscriptstyle M} }}=\{\boldsymbol{\theta}_{{ \mathcal{\scriptscriptstyle M} }}\in\mathbb{R}^m:|\boldsymbol{\theta}_{{ \mathcal{\scriptscriptstyle M} }}-\hat{\boldsymbol{\theta}}_{n,{ \mathcal{\scriptscriptstyle M} }}|_{\infty}\leq O_{{ \mathrm{p} }}(\nu)\}$, and $\tilde{\Lambda}_n(\boldsymbol{\theta}_{{ \mathcal{\scriptscriptstyle M} }})=\{\boldsymbol{\lambda}\in \mathbb{R}^m:\boldsymbol{\lambda}^{ \mathrm{\scriptscriptstyle \top} }{\mathbf f}_t^{{\mathbf A}_n}(\boldsymbol{\theta}_{{ \mathcal{\scriptscriptstyle M} }},\hat{\boldsymbol{\theta}}_{n,{ \mathcal{\scriptscriptstyle M} }^{\rm c}})\in \mathcal{V} \text{ for any }t\in[n]\}$. To justify this PPEL estimator, we impose one more assumption.
By definition, the population quantity ${\mathbf a}_k^0$ is the counterpart of the sample's ${\mathbf a}^n_k$. Condition (ref)(a) guarantees that these ${\mathbf a}_k^0$ are $L_1$-sparse and they are well estimated by ${\mathbf a}^n_k$ from ((ref)). Condition (ref)(b) ensures that ${\mathbf A}$ spans an $m$-dimensional linear space so the information provided by the moments concerning the $m$-dimensional parameters of interest does not degenerate.
Proposition (ref) provides the convergence rate of the PPEL estimator $\tilde{\boldsymbol{\theta}}_{{ \mathcal{\scriptscriptstyle M} }}$. If $L_n$ is a constant, $\tilde{\boldsymbol{\theta}}_{{ \mathcal{\scriptscriptstyle M} }}$ reaches the usual optimal rate of convergence $n^{-1/2}$. To find its asymptotic distribution, we must cope with the potential serial correlation in the projected moment functions ${\mathbf f}^{{\mathbf A}_n}(\cdot;\cdot)$. Define
as an estimator of the long-run covariance matrix, where
Here $\mathcal{K}(\cdot)$ is a symmetric kernel for the estimation of the long-run covariance matrix (newey1987simple, Andrews(1991)), inside of which lies the diverging bandwidth $h_n$. Condition (ref) consists of standard conditions on these kernels, which are satisfied by the widely used choices such as the Parzen kernel, the Tukey-Hanning kernel, and the QS kernel.
As the counterpart of (ref), we define
with $\widehat{{\mathbf V}}_{{\mathbf f}^{{\mathbf A}_n}}(\tilde{\boldsymbol{\theta}}_{{ \mathcal{\scriptscriptstyle M} }},\hat{\boldsymbol{\theta}}_{n,{ \mathcal{\scriptscriptstyle M} }^{\rm c}})=\mathbb{E}_n[\{{\mathbf f}^{{\mathbf A}_{n}}_t(\tilde{\boldsymbol{\theta}}_{{ \mathcal{\scriptscriptstyle M} }},\hat{\boldsymbol{\theta}}_{n,{ \mathcal{\scriptscriptstyle M} }^{\rm c}})\}^{\otimes2}]$. Now we are ready to state the asymptotic normality of the PPEL estimator.
Theorem (ref) establishes the asymptotic normality of $\tilde{\boldsymbol{\theta}}_{{ \mathcal{\scriptscriptstyle M} }}$. Unlike Theorem (ref), here the asymptotic distribution is well centered around zero. Standardized by a consistent estimate of the asymptotic variance, the limiting distribution is the desirable $N(0,1)$. In particular, when researchers are interested in the coefficient of one variable at a time, we can use ${\mathbf z}$ to pick out the coordinate and the left-hand side of the statement becomes a $t$-statistic; we then refer to a corresponding quantile of $N(0,1)$ to decide whether we shall reject the null hypothesis at a pre-specified test size.
Compared with Theorem (ref), Theorem (ref) further requires the tuning parameters $(\nu,\pi,\varsigma)$ and the bandwidth $h_n$ in the kernel to satisfy
and $\nu\varsigma\ll n^{-1/2}s^{-1/2}\ell_n^{-1/2}$ with $\ell_n \ll \min\{s^{-1}n^{-2/\gamma}\aleph_n^{-1},s^{-2}\aleph_n^{-1/2}, s^{-1/3}n^{-1/6}\aleph_n^{-2/3}, L_n^{-2} s^{-1/3} \aleph_n^{-2/3}\}$ and {$s\ll \min\{n^{\delta}, n^{3\varphi/(3\varphi+1)}L_n^{-6}\}$ for {$\delta= \min\{3\varphi/(6\varphi+2)-2/\gamma,3\varphi/(24\varphi+8)\}$}}. If we allow $s$ to grow in a polynomial order of $n$ in that $s\asymp n^{\kappa}$ satisfying $n^{\kappa}\ll n^{\psi}L_n^{-6} $, where $\kappa\in[0,{\delta})$, and $\psi=3\varphi/(3\varphi+1)$, then our proposed estimator $\tilde{\boldsymbol{\theta}}_{{ \mathcal{\scriptscriptstyle M} }}$ is asymptotically normal even if the number of moment conditions $r$ diverges at an exponential order in that $ \log r\ll \min \{ n^{\psi-4/\gamma-2\kappa}, n^{\psi-8\kappa},n^{(3\varphi-1)/(6\varphi+2)-\kappa}, n^{\psi-\kappa}L_n^{-6}, n^{\psi/3}L_n^{-\varphi}, n^{-\psi/(3\varphi)}\omega_n^{-2} \}, $ where $L_n\ll \min\{n^{1/(4\varphi)}(\log n)^{-1/\varphi},n^{\psi/6}, n^{1/7}\}$ and $\omega_n\ll \min\{L_n^{-3},n^{-1/(6\varphi+2)}\}$. The ranges of admissible $r$ and $p$ highlight the adaptivity to high-dimensional models.
To summarize the theoretical results, Theorem (ref) shows that the PEL estimator $\hat{\boldsymbol{\theta}}_n$ is a consistent estimator for the true parameter $\boldsymbol{\theta}_0$, and Theorem (ref) verifies its asymptotic normality of the nonzero components. However, the limiting normal distribution is not centered at zero due to the influence of the high dimensional moments, making it difficult to use for statistical inference. While most applied econometric use cases of inference focus on a low-dimensional parameter of interest, say $\boldsymbol{\theta}_{0,{ \mathcal{\scriptscriptstyle M} }}$, the PPEL estimator $\tilde{\boldsymbol{\theta}}_{{ \mathcal{\scriptscriptstyle M} }}$ projects out the influence of the nuisance parameters in $\mathcal{M}^{\mathrm c}$ from PEL, and it restores in Theorem (ref) the standard inferential procedure based on a zero-mean limiting normal distribution.
In this section we demonstrate our estimation and inference procedures in three important econometric models which fit seamlessly into our framework. The VAR model conventionally features a small cross section and a relatively long time dimension, thereby taking $n\to \infty$ only in the asymptotic framework. Such asymptotics fails to provide satisfactory approximation to the finite sample behavior when the cross section is non-trivial. The second example is the local projection method that is closely related to VAR, where the number of unknown parameters accumulates when the prediction horizon moves forward. The third example is the MGARCH that mimics VAR in modeling the dynamics and interaction of the volatility over a cross section.
Estimation for high-dimensional models involves meticulous tuning and optimization. PEL's $P_{1,\pi}(\cdot)$ and $P_{2,\nu}(\cdot)$ are chosen as the SCAD penalty and the $L_1$-norm (Lasso) penalty, respectively. We employ the interior-point method koh2007interior,koh2007efficient to efficiently solve the inner layer optimization for $\boldsymbol{\lambda}$ in (ref). For the outer layer, to handle the non-differentiability of the SCAD penalty at 0, we utilize a strategy that synthesizes adaptive moment estimation (ADAM) algorithm Kingma2014AdamAM with proximal gradient descent. The tuning parameters $\pi$ and $\nu$ are chosen by minimizing the BIC-type function
where $\hat{\boldsymbol{\theta}}_n^{(\nu,\pi)}$ is a local minimizer of (ref) for a given pair $(\nu,\pi)$, and $\text{df}(\hat{\boldsymbol{\theta}}_n^{(\nu,\pi)})$ and $\text{df}(\hat\boldsymbol{\lambda}(\hat{\boldsymbol{\theta}}_n^{(\nu,\pi)}))$ denote the number of nonzero elements in $\hat{\boldsymbol{\theta}}_n^{(\nu,\pi)}$ and $\hat\boldsymbol{\lambda}(\hat{\boldsymbol{\theta}}_n^{(\nu,\pi)})$, respectively. For the PPEL method, we solve (ref) for a properly chosen tuning parameter $\varsigma = 0.2n^{-1/3}$ which satisfies the theoretical assumptions. When we experiment with the Parzen kernel, the Tukey-Hanning kernel, and the QS kernel to compute the asymptotic variance of the PPEL estimator under the bandwidth $h_n=n^{1/5}$, we find the numerical results are robust to all kernels; we therefore report those under the Parzen kernel only. For each data generation process, we repeat $N=500$ times for computing $\text{MSE}=p^{-1}N^{-1}\sum_{i=1}^{N} |\hat\boldsymbol{\theta}^{(i)}-\boldsymbol{\theta}_0|_2^2$, $\text{Bias}^2 =p^{-1} |N^{-1}\sum_{i=1}^{N}\hat\boldsymbol{\theta}^{(i)}-\boldsymbol{\theta}_0|_2^2$ and $\text{Var} = \text{MSE} - \text{Bias}^2$, where $\hat\boldsymbol{\theta}^{(i)}$ is the estimate of $\boldsymbol{\theta}_0$ in the $i$-th repetition.
VAR is an off-the-shelf multivariate time series model. We first present a VAR($l$) in the following example to fix the notations.
{\bf Example 1: VAR($l$)}. A VAR of lag $l$ for a vector ${\mathbf z}_t =(z_{t,1},\ldots,z_{t,d_z})^{{ \mathrm{\scriptscriptstyle \top} }}\in\mathbb{R}^{d_z}$ follows
where ${\mathbf G}_1,\ldots,{\mathbf G}_l$ are $l$ coefficient matrices and $\boldsymbol{\varepsilon}_t$ is a white noise series. We collect ${\mathbf x}_t=({\mathbf z}_t^{{ \mathrm{\scriptscriptstyle \top} }},\ldots, {\mathbf z}_{t-l}^{{ \mathrm{\scriptscriptstyle \top} }})^{{ \mathrm{\scriptscriptstyle \top} }}$, $\boldsymbol{\theta}=\{{\rm vec}({\mathbf G}_1)^{{ \mathrm{\scriptscriptstyle \top} }},\ldots,{\rm vec}({\mathbf G}_l)^{{ \mathrm{\scriptscriptstyle \top} }}\}^{{ \mathrm{\scriptscriptstyle \top} }}$, and
Given $\boldsymbol{\varepsilon}_t$ is of zero mean and uncorrelated with ${\mathbf z}_{t-1},\ldots,{\mathbf z}_{t-l}$, we have $ \mathbb{E}\{{\mathbf g}({\mathbf x}_t;\boldsymbol{\theta}_0)\}={\bf0}$, which is over-identified for $\boldsymbol{\theta}_0$.
While AR(1) is the prototype of the univariate AR model, VAR(1) is the most used VAR specification in practice. We generate ${\mathbf z}_t$ from ${\mathbf z}_t = {\mathbf G}_1{\mathbf z}_{t-1}+\boldsymbol{\varepsilon}_t$ with $\boldsymbol{\varepsilon}_t\overset{{\rm i.i.d.}}{\sim}N({\mathbf 0}, \boldsymbol{\Sigma}_{\boldsymbol{\varepsilon}})$, where ${\mathbf G}_1$ is a $d_z\times d_z$ parameter matrix and $\boldsymbol{\Sigma}_{\boldsymbol{\varepsilon}}$ is specified as
The VAR(1) implies the moments $ \mathbb{E}(\boldsymbol{\varepsilon}_t)=\mathbb{E}({\mathbf z}_t - {\mathbf G}_1{\mathbf z}_{t-1})={\mathbf 0}\ $ and $\mathbb{E}(\boldsymbol{\varepsilon}_t\otimes{\mathbf z}_{t-1})=\mathbb{E}\{({\mathbf z}_t - {\mathbf G}_1{\mathbf z}_{t-1})\otimes{\mathbf z}_{t-1}\}={\mathbf 0}$.
We experiment with four combinations $(n,d_z) = (50,10), (80,10), (50,30)$ and $(80,30)$. The model complexity $d_z$ is comparable with the dimension in the empirical application, whereas $n$ here is sufficient to illustrate the performance. Under each $(n,d_z)$, we generate the coefficient matrix ${\mathbf G}_1$ with $10\%$ nonzero elements at random and rescale it to ensure that the process is stable with the signal-to-noise ratio $2:1$.
We compare the PEL estimator of $\boldsymbol{\theta}_0$ with the OLS estimator and basu2015regularized's $\ell_1$-LS estimator under their tuning procedure. We choose the OLS estimate as the initial value of PEL. The results are reported in Table (ref). Utilizing the sparsity of the parameters, the penalized methods significantly outperform OLS in various settings. When $d_z = 10$, PEL attains an edge over the competitive methods. When $d_z = 30$, the performances of the PEL estimator and the $\ell_1$-LS estimator are comparable, but PEL has the advantage of exhibiting lower bias.
As basu2015regularized does not provide an asymptotic distribution for inference, we can only focus on our PPEL method. The coverage frequency and median length of the PPEL-based CI for $\theta_{0,k}$ are tabulated in Table (ref). Without loss of generality, we simply choose $k$ as the index of the first nonzero element of $\boldsymbol{\theta}_0$. Three confidence levels, $90\%$, $95\%$ and $99\%$, are considered. Table (ref) shows that the coverage frequency approaches its nominal levels, and the median length of the CIs decreases as $n$ gets bigger with $d_z$ fixed. In addition, the median length of the CIs increases as $d_z$ gets bigger with $n$ fixed. These observations reflect the interaction between the complexity of the models and the sample sizes.
IRFs, interpreted as causal effects in dynamic modeling, are of central economic interest in macroeconomics. IRF can be elicited from the multi-equation VAR, or from a sequence of single-equation regressions. Let us continue with Example 1.
{\bf Example 2: Local Projection}. Without loss of generality, suppose that we are interested in the first variable $z_{t,1}$ of the vector ${\mathbf z}_t$. The VAR system's first equation is a one-period-ahead predictive regression $ z_{t,1} = {\mathbf G}_{1,\cdot} {\mathbf z}_{t-1} + \cdots + {\mathbf G}_{l,\cdot} {\mathbf z}_{t-l} + \varepsilon_{t,1}, $ where ${\mathbf G}_{k,\cdot}$ is the $k$-th row of ${\mathbf G}_k$. The key idea of LP is projecting the $h$-period-ahead $z_{t+h,1}$ to the past variables of time $(t-1, \ldots, t-l)$: $$ z_{t+h,1} = {\mathbf z}_{t-1}^{{ \mathrm{\scriptscriptstyle \top} }} \boldsymbol{\beta}_1^{(h)} + \cdots + {\mathbf z}_{t-l}^{{ \mathrm{\scriptscriptstyle \top} }} \boldsymbol{\beta}_l^{(h)} + \varepsilon^{(h)}_{t+h,1} $$ for each $h\in \{0\}\cup [H]$, where $H$ is the longest prediction horizon of interest, and $\{\boldsymbol{\beta}_1^{(h)}, \ldots, \boldsymbol{\beta}_l^{(h)} \}_{h=0}^H$ are functions of the VAR parameters by simple iterative substitution. It turns out that the sequence of $\{\boldsymbol{\beta}_1^{(h)}\}_{h=0}^H$ is the IRF Jorda2005. These regressions imply the $(H+1)\times l$ moment conditions
for $h\in \{0\}\cup [H]$ and $q\in [l]$, where $\boldsymbol{\theta} = \textrm{vec} ( \{\boldsymbol{\beta}_1^{(h)}, \ldots, \boldsymbol{\beta}_l^{(h)} \}_{h=0}^H )$.
Our simulation design mimics the empirical application based on RameyZubairy2018. They model two variables, $z_{t,1}$ for the GDP growth and $z_{t,2}$ for the growth of government spending, via a VAR(4) for “quarterly data”
We set the true parameters $ {\mathbf G}_1 = \bigg(
\bigg)$ and ${\mathbf b}_0 = (0.5, 0.5)^{{ \mathrm{\scriptscriptstyle \top} }}$, along with the sparse coefficients ${\mathbf b}_1={\mathbf 0}$, ${\mathbf G}_l={\mathbf 0}$ and ${\mathbf b}_l={\mathbf 0}$ for $l=2,3,4$, and the error term is generated from $ (\epsilon_{t,1}, \epsilon_{t,2})^{{ \mathrm{\scriptscriptstyle \top} }} \sim N({\mathbf 0}, \boldsymbol{\Sigma}_\epsilon)$ with $\boldsymbol{\Sigma}_\epsilon = \bigg(
\bigg)$. We use the real data of \cite{RameyZubairy2018}'s exogenous news shock for the variable ``$\mathrm{shock}_t$''. Without loss of generality, we look at the results for $z_{t,1}$. Given the simulated data, we estimate the linear regression
where ${\mathbf w}_{t-1}$ is a vector of control variables, and $\Psi_1^{(h)}(L)$ is a polynomial in the lag operator.
The LP model features the moment conditions as in (ref). The default implementation of LP is via OLS. We will compare the numerical performances of OLS and PEL/PPEL. Following the specification in RameyZubairy2018, we set $h=0,1,\ldots,20$ in (ref), the control variable vector ${\mathbf w}_{t-1}$ includes lags of $(z_{t,1}, z_{t,2}, \mathrm{shock}_t)$, and the number of the lags is $4$. Notice that here we have $p = 294$ parameters to handle, as each of the $h$-period ahead regressions involves an $\alpha^{(h)}_1$, a $\beta^{(h)}_1$, and $3 \times4 $ coefficients for lag terms, totaling $294 =21\times (1+1+12)$.
The simulation results are summarized in Table (ref) for $n=300$ and $500$, where we use “LP” to denote the standard LP's OLS estimator. Focusing on the IRF, here we report the coverage frequency and median CI length for each $\beta_1^{(h)}$. Under the panel “Estimation”, PEL enjoys much smaller variance and overall MSE than LP, showing the benefit of taking advantage of the sparsity. In terms of the coverage frequency for each coefficient $\beta_1^{(h)}$, here we report those of the even $h$ as the results under the odd $h$ have virtually no difference in terms of the patterns. The empirical coverage of LP is lower than the nominal counterpart, while PPEL (with long-run variance estimated under the Parzen kernel) performs much better. The observation that LP's coverage frequency does not improve as the sample size expands from $n=300$ to $n=500$ stems from actual data on military news shocks, as depicted in Figure (ref). Historically, a significant portion of America's major shocks --- 88% of the total 500 shocks by absolute magnitude --- happened within the initial 300 observations (prior to 1960). The standard deviation of these first 300 shocks is 0.077, whereas it notably drops to 0.012 for the following 200 shocks.
While VAR uses historical data to predict the means, MGARCH forecasts future volatility.
{\bf Example 3: MGARCH model}. Suppose a stochastic vector process ${\mathbf y}_t\in\mathbb{R}^{d_y}$ with $\mathbb{E}({\mathbf y}_t)={\bf0}$ follows EngleandKroner1995's MGARCH-BEKK(1,1) model
where ${\mathbf D}$ and ${\mathbf B}$ are $d_{y}\times d_{y}$ parameter matrices, ${\mathbf C}$ is a $d_{y}\times d_{y}$ triangular matrix. Denote $\mathcal{J}_{-\infty}^{t}$ as the $\sigma$-filed generated by $\{{\mathbf y}_s\}$ up to and including time $t$. Then by (ref), we have $ \mathbb{E}({\mathbf y}_{t}{\mathbf y}_{t}^{{ \mathrm{\scriptscriptstyle \top} }}|\,\mathcal{J}_{-\infty}^{t-2})={\mathbf C}^{{ \mathrm{\scriptscriptstyle \top} }}{\mathbf C}+{\mathbf D}\mathbb{E}({\mathbf y}_{t-1}{\mathbf y}_{t-1}^{{ \mathrm{\scriptscriptstyle \top} }}|\,\mathcal{J}_{-\infty}^{t-2}){\mathbf D}^{{ \mathrm{\scriptscriptstyle \top} }}+{\mathbf B}\mathbb{E}({\mathbf y}_{t-1}{\mathbf y}_{t-1}^{{ \mathrm{\scriptscriptstyle \top} }}|\,\mathcal{J}_{-\infty}^{t-2}){\mathbf B}^{{ \mathrm{\scriptscriptstyle \top} }} $, which implies
Let ${\mathbf q}^K({\mathbf y}_{t-2})=\{q^{K}_{1}({\mathbf y}_{t-2}),\ldots,q^{K}_{K}({\mathbf y}_{t-2})\}^{{ \mathrm{\scriptscriptstyle \top} }}$ denote a $K\times 1$ vector of known basis functions which, as $K\to \infty$, well approximate square integrable functions of ${\mathbf y}_{t-2}$, such as polynomial splines, B-splines, and power series. Then the conditional moment restrictions in (ref) lead to a large number of unconditional moments
where ${\mathbf x}_t = ({\mathbf y}_t^{{ \mathrm{\scriptscriptstyle \top} }},{\mathbf y}_{t-1}^{{ \mathrm{\scriptscriptstyle \top} }},{\mathbf y}_{t-2}^{{ \mathrm{\scriptscriptstyle \top} }})^{{ \mathrm{\scriptscriptstyle \top} }}$ and $\boldsymbol{\theta}_{0}=\{{\rm vech}({\mathbf C})^{{ \mathrm{\scriptscriptstyle \top} }},{\rm vec}({\mathbf D})^{{ \mathrm{\scriptscriptstyle \top} }}, {\rm vec}({\mathbf B})^{{ \mathrm{\scriptscriptstyle \top} }}\}^{{ \mathrm{\scriptscriptstyle \top} }}$. Moreover, (ref) also implies
Together, they produce the moment constraints ${\mathbf g}({\mathbf x}_t;\boldsymbol{\theta}_0)=\{{\mathbf g}_1({\mathbf x}_t;\boldsymbol{\theta}_0)^{{ \mathrm{\scriptscriptstyle \top} }}, {\mathbf g}_2({\mathbf x}_t;\boldsymbol{\theta}_0)^{{ \mathrm{\scriptscriptstyle \top} }}\}^{{ \mathrm{\scriptscriptstyle \top} }}$.
We experiment with MGARCH in (ref) under the following parameter matrices ${\mathbf C},{\mathbf D},{\mathbf B}$:
In our estimation, we specifically choose ${\mathbf q}^K({\mathbf y}_{t-2}) = (y_{t-2,1},\ldots, y_{t-2,5})^{{ \mathrm{\scriptscriptstyle \top} }}$ in (ref). For convenience in the simulation, we randomly choose a point near $\boldsymbol{\theta}_0$ as the initial value, $\boldsymbol{\theta}^0 = \boldsymbol{\theta}_0 + \boldsymbol{\varepsilon}$ where $\boldsymbol{\varepsilon} \sim N({\mathbf 0}, 0.5^2{\mathbf I}_p)$, for all three estimators: PEL, MLE and $\ell_1$-penalized MLE ($\ell_1$-MLE). Table (ref) summarizes the performances of the estimators. MLE and $\ell_1$-MLE fail to numerically converge in all settings --- in our experiments, MLE and $\ell_1$-MLE numerically break down even if the initial value is set as the true $\boldsymbol{\theta}_0$, for the dimension of the parameter here surpasses the capacity of these two methods. In sharp contrast, PEL maintains robustness. The MSE of PEL decreases as $d_y$ gets bigger under the same $n$, because the significantly higher parameter dimension $p = 2d_y^2 + d_y(d_y+1)/2$ of MGARCH enhances relative sparsity and thereby the effectiveness of the penalization. Moreover, Table (ref) shows that PPEL's coverage frequency aligns well with the nominal ones.
Given the reasonable performance of our PEL/PPEL procedure in the simulations, we apply this method to three real-data economic and financial applications.
Inflation is a key macroeconomic indicator which serves as a signal of broad economic conditions, helping policymakers, businesses, and households make informed decisions. Controlling inflation, as one of the Federal Reserve's dual mandates, is vital for maintaining economic stability and overall well-being. Inflation is measured by the growth rate of a price index. Besides the consumer price index (CPI), the personal consumption expenditures (PCE) price index is the Fed's primary measure for monetary policies. PCE has a wide coverage, with its sectoral indices breaking down overall price changes into sectors such as food, housing, energy, healthcare, and so on. Examining how prices evolve and co-move across sectors helps pinpoint the source of inflation.
Following the onset of the COVID-19 pandemic, the United States has seen considerable variations in its inflation rate. Notably, in December 2021, inflation peaked at 7.0%, marking the highest level in several decades, driven by surging demand and disruptions in the supply chain. This prompted the Federal Reserve to implement stricter monetary policies. By late 2023, the inflation rate had decreased to approximately 3.0%. The data utilized for this empirical analysis is sourced directly from the Bureau of Economic Analysis, which publishes the PCE. It includes quarter-to-quarter changes in 16 sectoral PCE indices that are seasonally adjusted, covering the period from the first quarter of 1959 to the third quarter of 2023 (259 quarters).
To understand the dynamics and sectoral spillover effects, we follow the implementation of our simulation in Section (ref) to fit a 16-sector VAR(1). The estimate of coefficient matrix is shown in Figure (ref). While some sectors have upstream-downstream relationships, other sectors are less connected. One salient feature is that most autoregressive coefficients (those on the diagonal) are non-zeros, which is a key driving force of the persistence of inflation. Secondly, the estimated matrix is overall quite sparse. Notice that V7: Gasoline and other energy goods is much more volatile than other sectors, due to weather conditions, geopolitical uncertainty, and occasional energy crises, and therefore the fitting mechanism delivers many more active coefficients along its row to reduce the magnitude of the corresponding residual.
We further follow DIEBOLD2014119's decomposition to measure the relative weight of the shocks. The $h$-step generalized variance decomposition matrix ${\mathbf D}^{g,h} = (d_{i,j}^{{g,h}})$ is $$ d_{i,j}^{g,h} = \frac{\sigma_{\varepsilon,j,j}^{-1}\sum_{\ell=0}^{h-1}(\boldsymbol{\iota}_i^{{ \mathrm{\scriptscriptstyle \top} }} {\mathbf G}_1^{\ell} \boldsymbol{\Sigma}_{\varepsilon}\boldsymbol{\iota}_j)^2 }{\sum_{\ell=0}^{h-1}\boldsymbol{\iota}_i^{{ \mathrm{\scriptscriptstyle \top} }}{\mathbf G}_1^{\ell}\boldsymbol{\Sigma}_{\varepsilon}{\mathbf G}_1^{\ell,{ \mathrm{\scriptscriptstyle \top} }}\boldsymbol{\iota}_i }\,, $$ where $\boldsymbol{\Sigma}_{\varepsilon}$ is the covariance matrix of the disturbance vector, $\sigma_{\varepsilon,j,j}$ is the $j$-th diagonal element of $\boldsymbol{\Sigma}_{\varepsilon}$, and $\boldsymbol{\iota}_i$ is the selection vector with one as the $i$-th element and zeros otherwise. We normalize each entry of the generalized variance decomposition matrix ${\mathbf D}^{g, h}$ by the row sum, denoted by $\widetilde{\mathbf D}^{g, h}$, and then calculate the column sums to measure the out-variance of each sector for different horizons $h$. The results are shown in Figure (ref). Each dot represents the relative fraction source of the variation that is originated from a particular sector over $h=1,\ldots, 10$. The two major sources of variation come from V7:Gasoline and other energy goods and V11:Transportation services. The former exhibits the highest variance among all sectors, as mentioned above, and the latter is tightly linked to the energy sector. They spread out the variations to other sectors via the interconnection of the network. This example illustrates that our PEL method allows us to work with the granular sectoral level indices to inspect the relative importance and dynamics, which is much richer than simply looking at an overall univariate inflation measure at the national level.
John Maynard Keynes introduced the government spending multiplier in his General Theory, which is foundational to his broader theories on fiscal policy and its role in managing economic activities, especially during recessions. Despite the importance of the concept, the measurement of the multiplier is by no means straightforward, because it is difficult to isolate it from confounding factors. An influential recent empirical study of the fiscal multiplier by RameyZubairy2018 uses quarterly data from 1889 to 2015. We follow the same specification to replicate their Figure 5, and then use our method to re-evaluate the IRFs. While our simulation in Section (ref) utilizes only the real data of news shock shown in Figure (ref), here in the empirical application we use the real time series of government spending and GDP. The LP specification in (ref) sets either government spending or GDP as the target variable, and $\text{shock}_t$ is again the military spending news. The control variable ${\mathbf w}_{t-1}$ includes four lags of the news shock, government spending, and GDP. A state-dependent alternative version of the model (ref) is also considered for the subsamples of high/low-unemployment states, respectively.
The empirical results are reported in Figure (ref). It consists of six subplots, each being the IRF of the fiscal multiplier of the news shock to government spending (the first row) and GDP (the second row). The three columns represent the full sample (the first column), the high unemployment periods (the second column, 181 observations), and the low-unemployment periods (the third column, 320 observations), following the original empirical study. The split of “bad times” and “good times” is to check the structural stability of the IRFs under distinctive states of the economy.
The two methods, RameyZubairy2018's LP estimated by OLS (which is the same as the original paper) and this paper's PPEL, return very close point estimates. The main message is that the fiscal multipliers are mostly much smaller than unity in American history through booms and recessions, so that counter-cyclical policies likely have dampened effects. However, there are visually salient gaps in terms of the interval estimates. Recall that our simulation in Section (ref) has shown that, despite the narrower length of CI, the coverage of LP may deviate from the nominal coverage frequency. It is echoed by Figure (ref) where PPEL has slightly wider CIs than LP for obtaining nominal CIs. In particular, in the second column of the subplots, the LP with a small $h =0, 1,\ldots,5$ has extremely narrow confidence intervals; given a sample size of 181 observations, it is surprising to see that the confidence intervals are even much narrower than those from the full sample of 501 as in the first column. The counterintuitive observation is likely the consequence of the low empirical coverage probability that deviates from the designated nominal coverage rate. On the other hand, the length of the confidence intervals from PPEL suitably reflects the statistical randomness that corresponds to the sample sizes. The evidence suggests that the latter provides more reasonable quantification of the underlying uncertainty.
In the last empirical application, we explore the volatility spillover effect of 16 stocks from China's banking industry, with data from the CSMAR database (\url{https://data.csmar.com/}). We investigate the daily log-returns of these stocks from January 2, 2024, to April 19, 2024. Table (ref) summarizes the descriptive statistics of the time series and their pairwise sample correlation. The sample means of all stocks are very close to 0, and the sample kurtosises of several stocks are much larger than 3. Moreover, all pairwise sample correlations are positive. These findings motivate us to employ the MGARCH-BEKK(1,1) model specified in (ref) to study the interconnection of volatility.
Recall that in the simulation in Section (ref) a disturbed true value was used as an initial value for PEL. Here for real data, we follow yao20241 to obtain an initial value. Specifically, we set the initial values for diagonal elements of ${\mathbf C},{\mathbf D}$ and ${\mathbf B}$ as square roots of estimated parameters of fitting each component series into a univariate GARCH(1,1) model, and the off-diagonal elements as random values sampled from $N(0, 0.5^2)$. Figure (ref) plots the estimated structure of $\widehat{\mathbf D}$ and $\widehat{\mathbf B}$. Each arrow corresponds to a nonzero estimated off-diagonal element of $\widehat{\mathbf D}$ or $\widehat{\mathbf B}$. Specifically, if $\widehat{\mathbf D}$'s $(i,j)$-th element $\hat{d}_{i,j}\neq 0$, the directional line shoots from $i$ to $j$.
To check the estimated patterns of connection, we further compute the numbers of outgoing links, incoming links, and net links of the estimated $\widehat{\mathbf D}$ graph, as suggested by dhaene2022volatility. Table (ref) shows that the major national banks are the sources of shocks. They include the four biggest ones (Industrial and Commercial Bank of China (GS), Agricultural Bank of China (NY), Bank of China (ZG), China Construction Bank (JS)), as well as Ping An Bank (PA), China Merchants Bank (ZS), and China CITIC Bank (ZX). These banks are central to China's financial system and play a significant role in sustaining the country's economic growth. On the other hand, several joint-stock commercial banks rank as top destinations of shocks. They mostly feature a higher level of market orientation and substantial exposure to market forces. These characteristics contribute to heightened sensitivity to market shocks and greater interconnectedness, resulting in more significant volatility spillover. The large national banks are the main focus in hedging the risk of China's banking industry. Our empirical exercise quantifies the volatility spillover effects, which can be of interest for both practitioners of risk management and policymakers.
This study investigates the PEL method in the analysis of multivariate time series models characterized by high-dimensional moments and parameters. These models are prevalent in the fields of economics and finance, as seen in VAR, local projection, and volatility models. We develop the marginal EL and demonstrate the consistency of PEL. Additionally, we introduce PPEL for the inference of low-dimensional parameters, which effectively removes the impact of nuisance parameters and maintains asymptotic normality centered at zero. Comprehensive Monte Carlo simulations provide evidence for the validity of our approach. We illustrate its practical application in three empirical examples using real data: the dynamics of the USA's inflation, the IRF of news shocks to government spending and GDP, and the stock price volatility spillovers among Chinese banks.
The PEL/PPEL method provides a flexible and adaptable strategy. It serves as a promising technique for regulating the high dimensionality in parameters and moments. This approach can be extended in various directions. For example, though this paper does not address time series models with endogeneity due to space constraints, instruments can be easily integrated into our framework using moment conditions.