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.
105,864 characters · 11 sections · 54 citation commands
Regularized Estimation of High-dimensional Factor-Augmented Vector Autoregressive (FAVAR) Models
{\bf Key words:} Model Identifiability; Compactness; Low-rank plus Sparse Decomposition; Finite-Sample Bounds\\ {\bf Source code:} \url{https://github.com/jhlinplus/High_dim_FAVAR_estimation}
There is a growing need in employing a large set of time series (variables) for modeling social or physical systems. For example, economic policy makers have concluded based on extensive empirical evidence sims1980macroeconomics,bernanke2005measuring,banbura2010large that large scale models of economic indicators provide improved forecasts, together with better estimates of how current economic shocks propagate into the future, which produces better guidance for policy actions. Another reason for considering large number of time series in social sciences is that key variables implied by theoretical models for policy decisions\footnote{such as the concept of output gap for monetary policy, the latter defined as the difference between the actual output of an economy and its potential output} are not directly observable, but related to a large number of other variables that collectively act as a good proxy of the unobservable key variables. In other domains such as genomics and neuroscience, advent of high throughput technologies have enabled researchers to obtain measurements on hundreds of genes from functional pathways of interest shojaie2010discovering or brain regions seth2015granger, thus allowing a more comprehensive modeling to gain insights into biological mechanisms of interest. There are two popular modeling paradigms for such large panel of time series, with the first being the Vector Autoregressive (VAR) model lutkepohl2005new and the second being the Dynamic Factor Model (DFM) stock2002forecasting,lutkepohl2014structural.
The VAR model has been the subject of extensive theoretical and empirical work primarily in econometrics, due to its relevance in macroeconomic and financial modeling. However, the number of model parameters increases quadratically with the number of time series included for each lag period considered, and this feature has limited its applicability since in many applications it is hard to obtain adequate number of time points for accurate estimation. Nevertheless, there is a recent body of technical work that leveraging {\em structured sparsity} and the corresponding regularized estimation framework has established results for consistent estimation of the VAR parameters under high dimensional scaling. basu2015estimation examined Lasso penalized Gaussian VAR models and proved consistency results, while at the same time providing technical tools useful for analysis of sparse models involving temporally dependent data. melnyk2016estimating extended the results to other regularizers, lin2017regularized to the inclusion of exogenous variables (the so-called VAR-X model in the econometrics literature), hall2016inference to models for count data and nicholson2017varx to the simultaneous estimation of time lags and model parameters. However, a key requirement for the theoretical developments is a spectral radius constraint that ensures the {\em stability} of the underlying VAR process basu2015estimation, lin2017regularized. For large VAR models, this constraint implies a smaller magnitude on average for all model parameters, which makes their estimation more challenging, unless one compensates with a higher level of sparsity. Nevertheless, very sparse VAR models may not be adequately informative, while their estimation requires larger penalties that in turn induce higher bias due to shrinkage, when the sample size stays fixed.
The DFM model aims to decompose a large number of time series into a few common latent factors and idiosyncratic components. The premise is that these common factors are the key drivers of the observed data, which themselves can exhibit temporal dynamics. They have been extensively used for forecasting purposes in economics stock2002forecasting, while their statistical properties have been studied in depth bai2008large. Despite their ability to handle very large number of time series, theoretically appealing properties and extensive use in empirical work in economics, DFMs aggregate the underlying time series and hence are not suitable for examining their individual cross-dependencies. Since in many applications researchers are primarily interested in understanding the interactions between key variables sims1980macroeconomics,stock2016dynamic, while accounting for the influence of many others so as to avoid model misspecification that leads to biased results, DFMs may not be the most appropriate model.
To that end, bernanke2005measuring proposed a “fusion" model, namely the Factor Augmented VAR, that aims to summarize the information contained in a large set of time series by a small number of factors and includes those in a standard VAR model. Specifically, let $\{F_t\}\in\mathbb{R}^{p_1}$ be the latent factor and $\{X_t\}\in\mathbb{R}^{p_2}$ the observed sets of variables, they jointly form a VAR system given by
In addition, there is a large panel of observed time series $Y_t\in\mathbb{R}^q$, whose current values are influenced by both $X_t$ and $F_t$; i.e., the calibration equation:
The primary variables of interest $X_t$ together with the unobserved factors $F_t$---both are assumed to have small and fixed dimensions---drives the dynamics of the system, and the factors are inferred from (ref).
Even in the low-dimensional setting ($p_2$ fixed), there is very limited theoretical work bai2016estimation on the FAVAR model and some work on identification restrictions for the model parameters bernanke2005measuring. However, the fixed dimensionality assumption is rather restrictive in many applications; in particular, the model has been extensively used in empirical work in economics and finance eickmeier2014understanding,caggiano2014uncertainty, yet customarily a very small size block $X_t$ is considered. For example, in bernanke2005measuring that introduces the FAVAR model, $X_t$ comprises of three “core" economic indicators (industrial production, consumer price index and the federal funds rate) and $Y_t$ of 120 other economic indicators. The VAR system is augmented by one factor summarizing the macroeconomic indicators, and the augmented system shows 7-lag time dependence that significantly increases the sample size requirement for estimation purposes. In a recent application, stock2016dynamic apply the FAVAR model to macroeconomics effects of oil supply shocks; the augmented VAR system consists of 8 times series (observed and latent), but due to the limitation in sample size to avoid non-stationarities ($T=120$) the lag of the model is fixed to 1. Hence, as argued in stock2016dynamic, there is growing need for large scale FAVAR models and this paper aims to examine their estimation and theoretical properties in high-dimensions, leveraging sparsity constraints on key model parameters.
The key contributions of this paper are twofold: (1) the introduction of an identifiability constraint compatible with the high-dimensional nature of the model, under sparsity assumptions on model parameters $\Gamma$ and $\{A^{(k)}\}$, and (2) the ensuing formulation of the optimization problem that leads to their estimators based on observational data and estimators' high-probability error bounds. At the technical level there are two sets of challenges that are successfully resolved: (i) the calibration equation involves both an observed set of covariates and a set of latent factors, and their interactions require careful handling to enable accurate estimation of the factors that constitute part of the input to the augmented VAR system and are crucial for estimating the transition matrix; and (ii) with the presence of a block of variables in the VAR system that are subject to error due to being estimated rather than directly observed, a number of new technical challenges emerge and they are compounded by the presence of temporal dependence. Note that for ease of presentation, the main technical developments are shown for Gaussian data (all noise processes in (ref) and (ref) are assumed to be Gaussian), but the key theoretical results are also established for sub-Gaussian and sub-exponential error processes; see Appendix C for a result of independent theoretical interest, even for the standard sparse VAR model.
\paragraph{Outline of the paper.} The remainder of the paper is organized as follows. In Section (ref), the model identifiability constraint is introduced, followed by formulation of the objective function to be optimized that obtains estimates of the model parameters. Theoretical properties of the proposed estimators, specifically, their high probability finite-sample error bounds, are investigated in Section (ref). Subsequently in Section (ref), we introduce an empirical implementation procedure for obtaining the estimates and present its performance evaluation based on synthetic data. An application of the model on interlinkages of commodity prices and the influence of world macroeconomic indicators on them is presented in Section (ref), while Section (ref) provides some concluding remarks. All proofs and other supplementary materials are deferred to Appendices.
\paragraph{Notations.} Throughout this paper, we use ${\vert\kern-0.25ex\vert\kern-0.25ex\vert A \vert\kern-0.25ex\vert\kern-0.25ex\vert}_{\cdot}$ to denote matrix norms for some generic matrix $A\in\mathbb{R}^{m\times n}$. For example, ${\vert\kern-0.25ex\vert\kern-0.25ex\vert A \vert\kern-0.25ex\vert\kern-0.25ex\vert}_1$ and ${\vert\kern-0.25ex\vert\kern-0.25ex\vert A \vert\kern-0.25ex\vert\kern-0.25ex\vert}_\infty$ respectively denote the matrix induced $1$-norm and infinity norm, ${\vert\kern-0.25ex\vert\kern-0.25ex\vert A \vert\kern-0.25ex\vert\kern-0.25ex\vert}_{\text{op}}$ the matrix operator norm and ${\vert\kern-0.25ex\vert\kern-0.25ex\vert A \vert\kern-0.25ex\vert\kern-0.25ex\vert}_{\text{F}}$ the Frobenius norm. Moreover, We use $\|A\|_1$ and $\|A\|_\infty$ respectively to denote the element-wise $1$-norm and infinity norm. For two matrices $A$ and $B$ of commensurate dimensions, denote their inner product by $\savebox{\@brx}{\(\m@th{\langle}\)} \mathopen{\copy\@brx\mkern2mu\kern-0.9\wd\@brx\usebox{\@brx}}A, B\savebox{\@brx}{\(\m@th{\rangle}\)} \mathclose{\copy\@brx\mkern2mu\kern-0.9\wd\@brx\usebox{\@brx}}= \text{tr}(A^\top B)$. Finally, we write $A\gtrsim B$ if there exists some absolute constant $c$ that is independent of the model parameters such that $A\geq cB$; and $A\asymp B$ if $A\gtrsim B$ and $B\gtrsim A$ hold simultaneously.
The FAVAR model proposed in bernanke2005measuring has the following two components, as seen in Section (ref): a system given in (ref) that describes the dynamics of the latent block $F_t\in\mathbb{R}^{p_1}$ and the observed block $X_t\in\mathbb{R}^{p_2}$ that jointly follow a stationary $\mathrm{VAR}(d)$ model (the “VAR equation"); and the model in (ref) that characterizes the contemporaneous dependence of the large observed informational series $Y_t\in\mathbb{R}^q$ as a linear function of $X_t$ and $F_t$ (the “calibration equation"). Further, $w^F_t$, $w^X_t$ and $e_t$ are all noise terms that are independent of the predictors, and we assume they are serially uncorrelated mean-zero Gaussian random vectors: $w^F_t\sim \mathcal{N}(0,\Sigma_w^F)$, $w^X_t\sim \mathcal{N}(0,\Sigma_w^X)$ and $e_t\sim \mathcal{N}(0,\Sigma_e)$. In this study we consider a potentially large VAR system that has many coordinates, hence in contrast to bernanke2005measuring and bai2016estimation where both $p_1$ and $p_2$ are fixed and small, we allow the size of the observed block, $p_2$, to be large\footnote{We do not impose the restriction that $p_2$ is smaller than the available sample size.} and to grow with the sample size; yet the size of the latent block, $p_1$, can not be too large and is still assumed fixed. Moreover, the size of the informational series, $q$, can also be large and grow with the sample size. Further, we assume that the transition matrices $\{A^{(i)}\}_{i=1}^d$ and the regression coefficient matrix $\Gamma$ are {\em sparse}. Finally, the factor loading matrix $\Lambda$ is assumed to be dense.
The latent nature of $F_t$ leads to the following observational equivalence across the following two models encoded by $(\Lambda,\Gamma)$ and $(\widetilde{\Lambda},\widetilde{\Gamma})$, respectively: for any invertible matrix $Q_1\in\mathbb{R}^{p_1\times p_1}$ and $Q_2\in\mathbb{R}^{p_1\times p_2}$,
where
In other words, the key model parameters $(\Lambda,\Gamma)$ and the latent factors $F_t$ are {\em not uniquely} identified, a known problem even in classical factor analysis anderson1958introduction. Thus, additional restrictions are required to overcome this indeterminacy, since there is an equivalence class parametrized by $(Q_1,Q_2)$ within which individual models are not mutually distinguishable based on observational data. For the FAVAR model, a total number of $p_1^2+p_1p_2$ restrictions are needed for unique identification of $\Lambda$, $\Gamma$ and $F_t$.
Various schemes have been proposed in the literature to address this issue. Specifically, bernanke2005measuring impose the necessary restrictions through the coefficient matrices of the calibration equation, requiring $\Lambda=\left[
\right]$ and $\Gamma_{[1:p_1],\cdot}=0$; that is, the upper $p_1\times p_1$ block of $\Lambda$ is set to the identity matrix and the first $p_1$ rows of $\Gamma$ to zero. \citet{bai2016estimation} consider different sets of restrictions that involve combinations of coefficients from the calibration equation and the noise term from the VAR equation. In the low-dimensional setting ($p_2$ fixed), one can proceed to estimate the parameters subject to these restrictions, by adopting either a single-step Bayesian likelihood approach \citep{bernanke2005measuring} or an orthogonal projection-based approach by profiling out $X_t$ \citep{bai2016estimation}. However, neither approach is applicable in high-dimensional settings, due to the growing dimension $p_2$ which would render a projection-based approach infeasible or add to the computational demands of a Bayesian procedure.
To overcome these issues in high-dimensional settings, we introduce an alternative identification scheme “IR$+$Compactness" that is compatible with the model specification and can also be seamlessly incorporated in the estimation procedure, leveraging sparsity of the regression coefficient $\Gamma$. Specifically, we first impose constraint (IR):
Note that (IR) imposes $p_1^2$ constraints but crucially not on the latent factors, given their subsequent utilization in the VAR system. Further, it yields uniquely identifiable $\Lambda$ and $F_t$, for any given product $\Lambda F_t$, and the indeterminacy incurred by $Q_1\in\mathbb{R}^{p_1\times p_1}$ in (ref) vanishes.
However, the issue is not fully resolved, since for any $Q_2\in\mathbb{R}^{p_1\times p_2}$, the following relationship holds:
where
All such models encoded by $(\check{F}_t,\check{\Gamma})$, form an equivalence class parametrized by $Q_2$ that specifies the transformation. We denote this equivalence class by $\mathcal{C}(Q_2)$. If $Q_2=O$, then $\mathcal{C}(Q_2)$ degenerates to a singleton that contains only the true data-generating model, which requires the imposition of $p_1p_2$ restrictions on primary model quantities. One applicable constraint out of theoretical consideration is to impose orthogonality on $X_t$ and $F_t$ --- it yields the necessary $p_1p_2$ restrictions; yet is excessively stringent and limits the appeal of the FAVAR model, while also being challenging to operationalize. Therefore as a good working alternative, we address the identifiability issue through a weaker constraint that effectively limits sufficiently the size of the $\mathcal{C}(Q_2)$.
To this end, let $\mathbf{X}\in\mathbb{R}^{n\times p_2}$, $\mathbf{Y}\in\mathbb{R}^{n\times q}$ and $\mathbf{F}\in\mathbb{R}^{n\times p_1}$ be centered data matrices whose rows are samples of $X_t$, $Y_t$ and the latent process $F_t$ respectively, and $\check{\mathbf{F}}$ is analogously defined. The characterization of $\mathcal{C}(Q_2)$ is through the sample versions of the underlying processes. Specifically, define the set of {\em factor hyperplanes} induced by $\mathcal{C}(Q_2)$ by
and we let $\Theta^\star$ denote the factor hyperplane associated with the true data-generating model, to distinguish it from some generic element in $\mathcal{S}(\check{\Theta})$ that is denoted by $\check{\Theta}$. Note that $\Theta^\star\in\mathcal{S}(\check{\Theta})$ and $\check{\Theta}$ coincides with $\Theta^\star$ when $Q_2=0$. Moreover, all elements in $\mathcal{S}(\check{\Theta})$ are at most of rank $p_1$, hence a low-rank component relative to their size $n\times q$. Next, in a similar spirit to negahban2012restricted, we define the following constrained set:
where $\varphi_{\mathcal{R}}(\Theta)$ is defined according to
and $\kappa(\mathcal{R}^*) :=\sup\nolimits_{\Theta\neq 0}\big( {\vert\kern-0.25ex\vert\kern-0.25ex\vert \Theta \vert\kern-0.25ex\vert\kern-0.25ex\vert}_{\text{F}}/\mathcal{R}^*(\Theta)\big)$ with $\mathcal{R}^*$ being the dual norm of some regularizer $\mathcal{R}$. Base on the above definition, $\varphi_{\mathcal{R}}(\Theta)$ captures the interaction between the factor space and the observed $\mathbf{X}$-space; the product $\kappa(\mathcal{R}^*)\mathcal{R}^*(\Theta)$ measures the spikiness of $\Theta$ w.r.t. $\mathcal{R}$, and in the case where $\mathcal{R}$ corresponds to the sparsity-induced $\ell_1$-norm which would be the setup of interest in this paper (see Section (ref)), $\mathcal{R}^*(\Theta) = \|\Theta\|_\infty$ and $\kappa(\mathcal{R}^*) = \sqrt{nq}$. With the definition of $\mathcal{S}_{\phi}(\check{\Theta})$, we impose the following compactness constraint on $\check{\Theta}$ to further encourage identifiability:
(Compactness) effectively limits the spikiness of all possible $\check{\Theta}$'s by imposing a {\em box constraint} through the dual norm corresponding to the sparsity regularizer, and for an arbitrary set of fixed realizations, it restricts the factor hyperplane set induced by $\mathcal{C}(Q_2)$ to its $\phi$-radius subset $\mathcal{S}_{\phi}(\check{\Theta})$. This in turn limits the size of the equivalence class $\mathcal{C}(Q_2)$ under consideration, since there is a one-to-one correspondence at the set level between $\mathcal{C}(Q_2)$ and the factor hyperplane set induced by it. This further implies that although the models encoded by $(F_t,\Gamma)$ and $(\check{F}_t,\check{\Gamma})$ may not be perfectly distinguishable based on observational data, at the population level the discordance between the two models can not be too large. It is worth pointing out that the bound $\phi(n,q)$ is allowed to grow, but at a much slower rate than the size of $\check{\Theta}$; specifically, we require $\phi(n,q)=o(\kappa(\mathcal{R}^*))$. For ease of presentation, we use $\phi$ to denote this bound henceforth and further note that it is in fact a constant in any finite sample setting.
In summary, our proposed identification scheme comprises of two parts: (IR) and (Compactness). The former provides exact identification within the factor hyperplane and narrows the scope of observationally equivalent models to $\mathcal{C}(Q_2)$, while the latter limits its size; and they jointly incur {\em approximate identification} of the true data generating model; and thus, for estimation purposes henceforth, it becomes adequate to focus on this restricted equivalence class, rather than its individual elements. The proposed scheme is suitable for the high-dimensional nature of the problem and can easily be incorporated in the formulation of the optimization problem for parameter estimation (see Section (ref)), which in turn yields estimates with tight error bounds (see Section (ref)).
Without loss of generality, we focus on the case where $d=1$ in subsequent technical developments, so that $Z_t:=(F^\top_t,X^\top_t)^\top$ follows a $\mathrm{VAR}(1)$ model $Z_t = AZ_{t-1} + W_t$:
The generalization to the $\mathrm{VAR}(d)$ $(d>1)$ case is straightforward since for any generic $\mathrm{VAR}(d)$ process satisfying $\mathcal{A}_d(L) Z_t = w_t$ where $\mathcal{A}_d(L):=\mathrm{I}-A^{(1)}L - \cdots - A^{(d)}L^d$, it can always be written in the form of a $\mathrm{VAR}(1)$ model for some $dp$-dimensional process $\widetilde{Z}_t$ lutkepohl2005new.
Based on the introduced model identification scheme (IR+Compactness), we propose the following procedure to estimate the FAVAR model, whose parameters include a sparse coefficient matrix $\Gamma$, a dense loading matrix $\Lambda$, and a sparse transition matrix $A$. Observed data matrices $\mathbf{X}$ and $\mathbf{Y}$ are identical to what have been previously defined, and to distinguish the responses from their lagged predictors when considering the VAR system, we let $\mathbf{X}_{n-1}:=[x_1,\dots,x_{n-1}]^\top$ denote the predictor matrix and $\mathbf{X}_{n}:=[x_2,\dots,x_n]^\top$ the response one; $\mathbf{F}_n, \mathbf{F}_{n-1}, \mathbf{Z}_n,\mathbf{Z}_{n-1}$ are analogously defined. Based on these notations, the sample versions of the VAR system and the calibration equation in (ref) and (ref) can be written as
We propose the following estimators obtained from a two-stage procedure for the coefficient matrices $\Lambda$, $\Gamma$ and subsequently the transition matrices $\{A_{ij}\}_{i,j=1,2}$.
In the presence of additional contemporaneous dependence amongst the coordinates for the error processes $w_t$, one may consider a maximum likelihood-based loss function, but the full estimation would require additional structural assumptions of $\Sigma_w$ (or its inverse) given the high dimensionality; we do not further elaborate in this study, since our prime interest is estimating the coefficient/transition matrices of the FAVAR model.
The formulation in (ref) based on the least squares loss function and the surrogate $\widehat{\mathbf{F}}$ is straightforward. However, the formulation for the calibration equation merits additional discussion. First, note that the factor hyperplane $\Theta$ has at most rank $p_1$ and therefore has low rank structure relative to its size $n\times q$. We impose a rank constraint in the estimation procedure to enforce such structure. Together with the (IR+Compactness) constraint introduced above, the objective then becomes to estimate accurately the parameters of a model within the equivalence class $\mathcal{C}(Q_2)$, in the sense that the estimate obtained by solving (ref) effectively corresponds to recovering an arbitrary $\check{\Theta},\check{\Theta}\in\mathcal{C}(Q_2)$; such an estimate, however, will be close to the true data generating $\Theta^\star$. Once this goal is achieved, this would enable accurate estimation of the transition matrix of the VAR system.
From an optimization perspective, the objective function admits a low-rank-plus-sparse decomposition and compactification is necessary for establishing statistical properties of the global optima in the absence of explicitly specifying the interaction structure between the low rank and the sparse blocks (or the spaces they live in). Note that the form of the compactness constraint is dictated by the statistical problem under consideration. For example, agarwal2012noisy study a multivariate regression problem, where the coefficient is decomposed to a sparse and a low rank block. In that setting, a compactness constraint is imposed through the entry-wise infinity norm bound of the low rank block. chandrasekaran2012latent study a graphical model with latent variables where the conditional concentration matrix is the parameter of interest. The marginal concentration matrix is decomposed to a sparse and a low rank block via the alignment of the Schur complement, and the compactness constraint is imposed on both blocks and manifests through the corresponding regularization terms in the resulting optimization problem. Hence, the compactness constraint takes different forms but ultimately serves the same goal, namely, to introduce an upper bound on the magnitude of the low rank--sparse block interaction, with the latter being an important component in analyzing the estimation errors. The compacteness constraint adopted for the FAVAR model serves a similar purpose, although the presence of temporal dependence introduces a number of additional technical challenges compared to the two aforementioned settings that consider independent and identically distributed data.
Finally, we remark that the model identification scheme (IR+Compactness) incorporated in the optimization problem as a constraint, enables us to establish high-probability error bounds (relative to the true data generating parameters/factors) for the proposed estimators, as shown next in Section (ref). Therefore, although (IR+Compactness) does not encompass the full $p_1^2+p_1p_2$ restrictions, it provides sufficient identifiability for estimation purposes.
In this section, we investigate the theoretical properties of the estimators proposed in Section (ref). We focus on formulations (ref) and (ref), whose global optima correspond to $(\widehat{\Theta},\widehat{\Gamma})$ and $\widehat{A}$, respectively.
Since (ref) relies not only on prime observable quantities (namely $X_t$), but also on estimated quantities from Stage I (namely $\widehat{\mathbf{F}}$), the analysis requires a careful examination of how the estimation error in the factor propagates to that of $\widehat{A}$. We start by outlining a road map of our proof strategy together with a number of regularity conditions needed in subsequent developments. Section (ref) establishes error bounds for $\widehat{\Gamma}$, $\widehat{\Theta}$ \footnote{Consequently, the error bounds of $\widehat{\mathbf{F}}$ and $\widehat{\Lambda}$ under (IR) are also obtained.}and $\widehat{A}$ under certain regularity conditions and employing suitable choices of the tuning parameters, for {\em deterministic realizations} from the underlying observable processes. Specifically when considering the error bound of $\widehat{A}$, the error of the plug-in estimate $\widehat{\mathbf{F}}$ is assumed non-random and given. Subsequently, Section (ref) examines the probability of the events in which the regularity conditions are satisfied for {\em random realizations}, and further establishes high-probability upper bounds for quantities to which the tuning parameters need to conform. Finally, the high-probability finite sample error bounds for the estimates obtained based on random realizations of the data generating processes readily follow after properly aligning the conditioning arguments, and the results are presented in Section (ref). All proofs are deferred to Appendices (ref) and (ref).
\paragraph{Additional notations.} Throughout, we use superscript $\star$ to denote the true value of the parameters of interest, and $\Delta$ for errors of the estimators; e.g., $\Delta_{A}=\widehat{A}-A^\star$. For sample quantities (e.g., $\mathbf{X}$ and $\mathbf{F}$) and their corresponding error (e.g., $\Delta_{\mathbf{F}}$), we use subscript $(n-1)$ to denote their first $n-1$ rows. We let $S_{\mathbf{E}}:=\tfrac{1}{n}\mathbf{E}^\top\mathbf{E}$ denote the sample covariance matrix of $\mathbf{E}$ and the sample covariance of other quantities are analogously defined. Additionally, denote the density level of $\Gamma^\star$ by $s_{\Gamma^\star}:=\|\Gamma^\star\|_0$, and that of $A^\star$ by $s_{A^\star}$.
\paragraph{A road map for establishing consistency results.} As previously mentioned, the key steps are:
In Part 1, note that the first-stage estimators obtained from the calibration equation are based on observed data and thus the regularity conditions needed are imposed on (functions of) the observed samples. On the other hand, the second-stage estimator relies on the plugged-in first-stage estimates that have bounded errors; therefore, the analysis is carried out in an analogous manner to problems involving error-in-variables. Specifically, the required regularity conditions on quantities appearing in the optimization (ref) involve the error of the first stage estimates, with the latter assumed fixed. In Part 2, the focus shifts to the probability of the regularity conditions being satisfied under random realizations, again starting from the first stage estimates, with the aid of Gaussian concentration inequalities and proper accounting for temporal dependence. Once the required regularity conditions are shown to hold with high probability, combining the results established in Part 1 for deterministic realizations, the high-probability error bounds for $\widehat{\Theta}$ and $\widehat{\Gamma}$ are established. The high-probability error bound of the estimated factors readily follows, which ensures that the variables which Stage II estimates rely upon are sufficiently accurate with high probability. Based on the latter result, the regularity conditions required for the Stage II estimates are then verified to hold with high probability at a certain rate. In the FAVAR model, since the estimation of the VAR equation is based on quantities among which one block is subject to error, to obtain an accurate estimate of the transition matrix requires more stringent conditions on population quantities (e.g., extremes of the spectrum), so that the regularity conditions hold with high probability. In essence, the joint process $Z_t$ need to be adequately “regular" in order to get good estimates of the transition matrix , vis-a-vis the case of the standard VAR model where all variables are directly observed.
Next, we introduce the following key concepts that are widely used in establishing theoretical properties of high-dimensional regularized $M$-estimators negahban2012unified,loh2012high, as well as quantities that are related to processes exhibiting temporal dependence basu2015estimation.
Under the current model setup, however, the exact form of the deviation bound becomes more involved and requires proper modifications to incorporate quantities associated with the factor hyperplane, as seen in Proposition (ref).
\medbreak We start by providing error bounds for $\widehat{\Gamma}$ and $\widehat{\Theta}$, as well as those of the corresponding $\widehat{\mathbf{F}}$ and $\widehat{\Lambda}$ extracted under (IR). For the optimization problem given in (ref), we assume that $r\geq p_1$ and $\phi$ is always compatible with the true data generating mechanism, so that $\Theta^\star$ is always feasible. To this end, the error bounds of $\widehat{\Theta}$ and $\widehat{\Gamma}$ for deterministic realizations crucially rely on two components: (i) $\mathbf{X}$ satisfying the RSC condition with curvature $\alpha_{\text{RSC}}^{\mathbf{X}}$; and (ii) the tuning parameter $\lambda_\Gamma$ being chosen in accordance with the deviation bound condition that is associated with the interaction between $\mathbf{X}$ and $\mathbf{E}$, the strength of the noise, and the interaction between the space spanned by the factor hyperplane and the observed $\mathbf{X}$. Upon the satisfaction of these conditions, the error bounds of $\widehat{\Theta}$ and $\widehat{\Gamma}$ are given by
and these conditions hold with high probability for random realizations of $X_t$ and $Y_t$. Since $\widehat{\mathbf{F}}$ is the first $p_1$ columns of $\widehat{\Theta}$, it possesses an error bound of the similar form.
Next, we briefly sketch the error bounds of $\widehat{A}$. For the optimization in (ref), for deterministic realizations, the results in basu2015estimation can be applied with the corresponding RSC condition and deviation condition imposed on quantities associated with $\widehat{\mathbf{Z}}_n$ and $\widehat{\mathbf{Z}}_{n-1}$, and the error for $\widehat{A}$ is in the form of
Then, for random realizations, assuming $\Delta_{\mathbf{F}}$ known and non-random, to satisfy the corresponding regularity conditions, we additionally require that the following functional involving the spectral density of the underlying joint process $Z_t$ exhibits adequate curvature, that is, $\mathfrak{m}(f_Z)/\sqrt{\mathcal{M}(f_Z)} > c_0 h_1(\Delta_{\mathbf{F}_{n-1}})$ for constant $c_0$ and some function $h_1$ of the error $\Delta_{\mathbf{F}_{n-1}}$ that captures its magnitude. Moreover, the deviation bound is of the form $h_2(\Delta_{\mathbf{F}})$, which can be viewed as another function of the error\footnote{note the deviation bound in principle also depends on other population quantities such as $\mathfrak{m}(f_Z)$, $\mathcal{M}(f_Z)$, $\Lambda_{\max}(\Sigma_w)$ etc.}. Further, since $\Delta_{\mathbf{F}}$ is bounded with high probability from the analysis in Stage I, it will be established that $h_1(\Delta_{\mathbf{F}})$ and $h_2(\Delta_{\mathbf{F}})$ are both upper bounded at a certain rate, thus ensuring that the RSC condition and the deviation conditions can both be satisfied unconditionally, by properly choosing the required constants.
Proposition (ref) below gives the error bounds for the estimators in (ref), assuming certain regularity conditions hold for deterministic realizations of the processes $X_t$ and $Y_t$, upon suitable choice of the regularization parameters.
Based on Proposition (ref), under fixed realizations of $X_t$ and $Y_t$, the error bounds of $\widehat{\Gamma}$ and $\widehat{\Theta}$ are established. Using these Stage I estimates and the IR condition, estimates of the factors and their loadings can be calculated. In particular, since $\Delta_{\mathbf{F}}$ corresponds to the first $p_1$ columns of $\Delta_{\Theta}$, the above bound automatically holds for $\Delta_{\mathbf{F}}$. Further, the following lemma provides the relative error of the estimated $\Lambda$ under (IR) and the condition on $\Lambda^{1/2}_{\max}(S_{\mathbf{F}})$, with the latter translating to the requirement that the leading signal of $\mathbf{F}$ overrules the averaged row error of $\Delta_{\Theta}$.
Up to this point, error bounds have been obtained for all the parameters in the calibration equation. The following proposition establishes the error bound for the estimator obtained from solving (ref), based on observed $\mathbf{X}$ and estimated $\widehat{\mathbf{F}}$, and assuming $\Delta_{\mathbf{F}}$ is fixed.
Note that Proposition (ref) applies the results in basu2015estimation to the setting in this study, where Stage II estimation of the transition matrix is based on $\widehat{\mathbf{Z}}_n$ and $\widehat{\mathbf{Z}}_{n-1}$; consequently, the regularity conditions should be imposed on corresponding quantities associated with $\widehat{\mathbf{Z}}_{n}$ and $\widehat{\mathbf{Z}}_{n-1}$.
Propositions (ref) and (ref) give finite sample error bounds for the estimators of the parameters obtained by solving optimization problems (ref) and (ref) based on fixed realizations of the observable processes $X_t$ and $Y_t$, and the regularity conditions outlined. Next, we examine and verify these conditions for random realizations of the processes, to establish high probability error bounds for these estimators.
We provide high probability bounds or concentrations for the quantities associated with the required regularity conditions, for random realizations of $X_t$ and $Y_t$. Specifically, we note that when $X_t$ is considered separately from the joint system, it follows a high-dimensional VAR-X model lin2017regularized
whose spectrum $f_X(\omega)$ satisfies
where $\mathcal{A}_{X}(L):=\mathrm{I}-A_{22}L$. Similar properties hold for $F_t$. Throughout, we assume $\{X_t\},\{F_t\}$ and $\{Y_t\}$ are all mean-zero stable Gaussian processes.
Lemmas (ref) to (ref) respectively verify the RSC condition associated with $\mathbf{X}$ and establish the high probability bounds for $\|\mathbf{X}^\top\mathbf{E}/n\|_\infty$, $\Lambda_{\max}(S_{\mathbf{E}})$ and $\Lambda_{\max}(S_{\mathbf{X}})$.
In the next two lemmas, we verify the RSC condition for random realizations of $\widehat{\mathbf{Z}}_{n-1}$ and obtain the high probability bound $C(n,p_1,p_2)$ for $\|\widehat{\mathbf{Z}}_{n-1}^\top\big( \widehat{\mathbf{Z}}_n - \widehat{\mathbf{Z}}_{n-1}(A^\star)^\top \big)/n\|_\infty$, with the underlying truth $\mathbf{F}$ being random but the error $\Delta_{\mathbf{F}}$ non-random. Note that this can be equivalently viewed as a {\em conditional} RSC condition and deviation bound, when conditioning on some fixed $\Delta_{\mathbf{F}}$.
Given the results in Sections (ref) and (ref), we provide next high probability error bounds for the estimates, obtained by solving the optimization problems in (ref) and (ref) based on random snapshots from the underlying processes $X_t$ and $Y_t$.
Theorem (ref) combines the results in Proposition (ref) and Lemmas (ref) to (ref) and provides the high probability error bound of the estimates, when $\widehat{\Theta}$ and $\widehat{\Gamma}$ are estimated based on random realizations from the observable processes $X_t$ and $Y_t$, with the latter driven by both $X_t$ and the latent $F_t$.
Note that the above bound also holds if we replace $\Delta_{\Theta}$ by $\Delta_{\mathbf{F}}$ under (IR). Next, using the results in Proposition (ref), Lemmas (ref) and (ref) and combine the bound in Theorem (ref), we establish a high probability error bound for the estimated $\widehat{A}$ in Theorem (ref).
We first discuss implementation issues of the proposed problem formulation for the high-dimensional FAVAR model. Specifically, the formulation requires imposing the compactness constraint for identifiability purposes and for obtaining the necessary statistical guarantees for the estimates of the model parameters. However, the value $\phi$ in the compactness constraint is hard to calibrate in any real data set. Hence, in the implementation we relax this constraint and assess the performance of the algorithm. Due to its importance in constraining the size of the equivalence class $\mathcal{C}(Q_2)$, we examine in Appendix (ref) certain relatively extreme settings where the proposed relaxation fails to provide accurate estimates of the model parameters.
\paragraph{Implementation.} The following relaxation of (ref) is used in practice:
which leads to Algorithm (ref).
The implementation of Stage I requires the pair of tuning parameters $(\lambda_\Gamma,r)$ as input, and the choice of $r$ is particularly critical since it determines the effective size of the latent block. In our implementation, we select the optimal pair based on the Panel Information Criterion (PIC) proposed in ando2015selecting, which searches for $(\lambda_\Gamma,r)$ over a lattice that minimizes
where $\widehat{\sigma}^2 = \tfrac{1}{nq}{\vert\kern-0.25ex\vert\kern-0.25ex\vert \mathbf{Y}-\widehat{\Theta} - \mathbf{X}\widehat{\Gamma}^\top \vert\kern-0.25ex\vert\kern-0.25ex\vert}_{\text{F}}^2$. Analogously, the implementation of Stage II requires $\lambda_A$ as input, and we select $\lambda_A$ over a grid of values that minimizes the Bayesian Information Criterion (BIC):
where $\text{RSS}_i:= \|(\mathbf{X}_{n})_{\cdot i}-\mathbf{X}_{n-1}\widehat{A}^\top_{i\cdot}\|^2$ is the residual sum of square of the $i^{\text{th}}$ regression. Extensive numerical work shows that these two criteria select very satisfactory values for the tuning parameters, which in turn yield highly accurate estimates of the model parameters.
\paragraph{Simulation setup.} Throughout, we assume $\Sigma_w^X$, $\Sigma_X^F$ and $\Sigma_e$ are all diagonal matrices, and the sample size is fixed at 200, unless otherwise specified. We first generate samples of $F_t\in\mathbb{R}^{p_1}$ and $X_t\in\mathbb{R}^{p_2}$ recursively according to the $\mathrm{VAR}(d)$ model in (ref), and then the samples of $Y_t\in\mathbb{R}^q$ are generated according to the linear model given in (ref). Specifically, (IR) is imposed on the true value of the parameter, hence $\Lambda^\star$ that is used for generating $Y_t$ always satisfies the restriction $\Lambda=\left[
\right]$. Unless otherwise specified, all error terms are generated according to some mean-zero Gaussian distribution.
For the calibration equation, the density level of the sparse coefficient matrix $\Gamma\in\mathbb{R}^{q\times p_2}$ is fixed at $5/p_2$ for each regression; thus, each $Y_t$ coordinate is affected by 5 series (coordinates) from the $X_t$ block on average. The bottom $(q-p_1)\times p_1$ block of the loading matrix $\Lambda\in\mathbb{R}^{q\times p_1}$ is dense. The magnitude of nonzero entries of $\Gamma$ and that of entries of $\Lambda$ may vary to capture different levels of signal contributions to $Y_t$, and we adjust the standard deviation of $e_t$ to maintain the desired level of the signal-to-noise ratio for $Y_t$ (averaged across all coordinates).
For the transition matrix $A$ of the VAR equation, the density for each of its component block $\{A_{ij}\}_{i,j=1,2}$ varies across settings, so as to capture different levels of the influence from the lagged value of the latent block $F_{t}$ on the observed $X_t$. Note that to ensure stability of the VAR system, the spectral radius of $A$, $\varrho(A)$, needs to be smaller than 1. In particular, when a $\mathrm{VAR}(d)~(d>1)$ system is considered, we need to ensure that the spectral radius of $\widetilde{A}$ is smaller than 1\footnote{In practice, this can be achieved by first generating $A^{(1)},\dots,A^{(d)}$, align them in $\widetilde{A}_{\text{initial}}$ and obtain the scale factor $\zeta:=\varrho_{\text{target}}/\varrho(\widetilde{A}_{\text{initial}})$, then scale $A^{(i)}$ by $\zeta^i$. The validity of this procedure follows from simple algebraic manipulations.}, where we let $p=p_1+p_2$ and
Table (ref) lists the simulation settings and their parameter setup.
Specifically, in settings A1\,--\,A4, $(F^\top_t,X^\top_t)^\top$ jointly follows a $\mathrm{VAR}(1)$ model. The (average) signal-to-noise ratio for each regression of $Y_t$ is 1.5. For settings A1 and A2, the transition matrix $A$ is uniformly sparse, with A2 corresponding to a larger system; for settings A3 and A4, we increase the density level (the proportion of nonzero entries) for the transition matrices that govern the effect of $F_{t-1}$ on $F_{t}$ and $X_{t}$. In particular, for setting A4, we consider a large system with 500 coordinates in $X_t$, and the factor effect is almost pervasive on these coordinates (through the lags), as the density level of $A_{21}$ is set at 0.8. Settings B1, B2 and B3 consider settings with more lags ($d=2$ and $d=4$, respectively), and to compensate for the higher level of correlation between $F_t$ and $X_t$, we elevate the signal-to-noise for each regression of $Y_t$ to 2. For B1, the transition matrices for both lags ($A^{(1)}$ and $A^{(2)}$) have uniform sparsity patterns, with $A^{(2)}$ being slightly more sparse compared to $A^{(1)}$; for B2, the transition matrices for the first two lags have higher density in the component that governs the $F_{t-i}\rightarrow X_t$ cross effect, and those for the last two lags have uniform sparsity. B3 has approximately the same scale as observed in real data, and due to a small $p_2$, the system exhibits a higher sparsity level in general. In settings C1\,--\,C4, the error terms of the VAR system are generated from distributions with tails heavier than a Gaussian (e.g. $t$-distributions, squares of Gaussian which have sub-exponential tails), and the joint process $(F_t',X_t)'$ will be heavy-tailed as a result of the recursive data generating mechanism.
\paragraph{Performance evaluation.} We consider both the estimation and the forecasting performance of the proposed estimation procedure. The performance metrics used for estimation are sensitivity (SEN), specificity (SPC) and the relative error in Frobenius norm (Err) for the sparse components (transition matrices $A$ and the coefficient matrix $\Gamma$), defined as
We also track the estimated size of the latent component (i.e., the rank constraint in (ref), jointly with $\lambda_{\Gamma}$ is selected by PIC), as well as the relative errors of $\widehat{\Theta}$, $\widehat{\mathbf{F}}$ and $\widehat{\Lambda}$. For forecasting, we focus on evaluating the $h$-step-ahead predictions for the $X_t$ block. Specifically, for settings where the VAR system is 1-lag dependent (A1--A4, C1), we consider $h=1$; for settings where the VAR system has more lag dependencies (B1--B3, C2--C4), we consider $h=1,2$. We use the same benchmark model as in banbura2010large which is based on a special case of the Minnesota prior distribution litterman1986forecasting, so that the for any generic time series $X_t\in\mathbb{R}^p$, each of its coordinates $j=1,\dots,p$ follows a centered random walk:
For each forecast $\widehat{x}_{T+h}$, its performance is evaluated based on the following two measures:
where $\text{rel-err}$ measures the $\ell_2$ norm of the relative error of the forecast to the true value; whereas for $\text{rel-err-ratio}$, it measures the ratio between the relative error of the forecast and the above described benchmark. In particular, its numerator and denominator respectively capture the averaged relative error of all coordinates of the forecast $\widehat{x}_{T+h}$ and that of the benchmark $\widetilde{x}_{T+h}$ that evolves according to (ref), while the ratio measures how much the forecast based on the proposed FAVAR model outperforms $(<1)$ or under-performs $(>1)$ compared to the benchmark.
All tabulated results are based on the average of 50 replications. Tables (ref), (ref) and (ref), respectively, depict the performance of the estimates of the parameters in the calibration and the VAR equations, as well as the forecasting performance under the settings considered.
Based on the results listed in Tables (ref) and (ref), we notice that in all settings, the parameters in the calibration equation $\widehat{\Theta}$ and $\widehat{\Gamma}$ are well estimated, while the rank slightly underestimated. Further, the SEN and SPC measures of $\widehat{\Gamma}$ show excellent performance regarding support recovery. It is worth pointing out that the estimation accuracy of the parameters in the calibration equation strongly depends on the signal-to-noise ratio of $Y_t$. In particular, if the signal-to-noise ratio in A1-A4 is increased to 1.8, the rank is always correctly selected by PIC, and the estimation relative error of $\widehat{\Theta}$ further decreases(results omitted for space considerations)\footnote{This also comes up when comparing the relative error of $\widehat{\Theta}$ in the A1-A4 settings to that in the B1-B2 ones, where the latter two have a higher SNR.}. Under the given IR, we decompose the estimated factor hyperplane into the factor block and its loadings. The results show that both quantities exhibit a higher relative error compared to that of the factor hyperplane. Of note, the loadings estimates exhibit a lot of variability as indicated by the high standard deviation in the Table.
Regarding the estimates in the VAR equation, for settings A1, A2 and B1 that are characterized by an adequate degree of sparsity, the recovery of the skeleton of the transition matrices is very good. However, performance deteriorates if the latent factor becomes “more pervasive" (settings A3 and A4), which translates to the $A_{21}$ block having lower sparsity. On the other hand, this does not have much impact on the recovery of the $A_{22}$ sub-block, as for these two settings, SEN and SPC of $A_{22}$ still remain at a high level. For settings with more lags, performance deteriorates (as expected) although SEN and SPC remain fairly satisfactory. On the other hand, the relative error of the transition matrices increases markedly. Nevertheless, the estimates of the first lag transition matrix is better than the remaining ones. Further, the results indicate that smaller size VAR systems (B3) exhibit better performance than larger ones. Finally, in terms of forecasting (results depicted in Table (ref)), the one-step-ahead forecasting value yields approximately 50% to 90% rel-err (compared to the truth), depending on the specific setting and the actual SNR, while it outperforms the forecast of the benchmark by around 40% (based on the $\text{rel-err-ratio}$ measure). Of note, the $2$-step-ahead forecasting value for settings with more lags outperforms the benchmark by an even wider margin with the rel-err-ratio decreasing to less than 0.3.
Finally, the proposed methodology is robust in the presence of heavier than Gaussian tails in the VAR processes. Further, note that in setting C3 wherein the temporal dependence is strong and the error terms are generated according to a sub-exponential distribution, the performance of the estimated transition matrices deteriorates significantly, as expected from the theoretical results outlined in Appendix C. Nevertheless, with proper compensation in terms of sample size (setting C4), the performance improves markedly.
Interlinkages between commodity prices represent an active research area in economics and have been a source of concern for policymakers. Commodity prices, unlike stocks and bonds, are determined more strongly by global demand and supply considerations. Nevertheless, other factors are also at play as outlined next. The key ones are: (i) the state of the global macro-economy and the state of the business cycle that manifest themselves as direct demand for commodities; (ii) monetary policy, specifically, interest rates that impact the opportunity cost for holding inventories, as well as having an impact on investment and hence production capacity that subsequently contribute to changes in supply and demand in the market; and (iii) the relative performance of other asset classes through portfolio allocation frankel2008effect,frankel2014effects. We employ the FAVAR model and the proposed estimation method to investigate interlinkages amongst major commodity prices. The $X_t$ block corresponds to the set of commodity prices of interest, while the $Y_t$ block contains representative indicators for the global economic environment. We extract the factors $F_t$ based on the calibration equation and then consider the augmented VAR system of $(F_t,X_t)$, so that the estimated interlinkages amongst commodity prices are based on a larger information set that takes into account broader economic activities.
\paragraph{Data.} The commodity price data ($X_t$) are retrieved from the International Monetary Fund, comprising of 16 commodity prices in the following categories: Metal, Energy (oil) and Agricultural. The set of economic indicators ($Y_t$) contain core macroeconomic variables and stock market composite indices from major economic entities including China, EU, Japan, UK and US, with a total number of 54 indicators. Specifically, the macroeconomic variables primarily account for: Output & Income (e.g. industrial production index), Labor Market (unemployment), Money & Credit (e.g. M2), Interest & Exchange Rate (e.g. Fed Funds Rate and the effective exchange rate), and Price Index (e.g. CPI). For variables that reflect interest rates, we use both the short-term interest rate such as 6-month LIBOR, and the 10-year T-bond yields from the secondary market. Further, to ensure stationarity of the time series, we take the difference of the logarithm for $X_t$; for $Y_t$, we apply the same transformation as proposed in stock2002forecasting. A complete list of the commodity prices and economic indicators used in this study is provided in Appendix (ref). For all time series considered, we use monthly data spanning the January 2001 to December 2016 period. Further, based on previous empirical findings in the literature related to the global financial crisis of 2008 stock2017twenty, we break the analysis into the following three sub-periods stock2017twenty: pre-crisis (2001--2006), crisis (2007--2010) and post-crisis (2011--2016), each having sample size (available time points) 72, 48, and 72, respectively\footnote{For each individual time series, we test for its normality using data spanning the pre-crisis, crisis, and post-crisis periods, respectively. Based on the Shapiro-Wilk test, the null hypothesis of normality is not rejected for selected time series (e.g., ALUMINUM) and rejected for others (e.g., OIL). However, when testing for multivariate normality of the joint distribution of all time series resp. across the three periods, we fail to reject the null hypothesis. The latter result may be due to inadequate power of the test given the relatively small sample size.}.
We apply the same estimation procedure for each of the above three sub-periods. Starting with the calibration equation, we estimate the factor hyperplane $\Theta$ and the sparse regression coefficient matrix $\Gamma$, then extract the factors based on the estimated factor hyperplane under the (IR) condition. For each of the three sub-periods, 4, 3, and 3 factors are respectively identified based on the PIC criterion, with the key variable loadings (collapsed into categories) on each extracted factor listed in Table (ref), after adjusting for $\Gamma X_t$.
Based on the composition of the factors, we note that the factors summarize both the macroeconomic environment and also capture information from the secondary market (bond & equity return), as suggested by economic analysis of potential contributors to commodity price movements frankel2008effect,frankel2014effects. Hence, the obtained factors summarize the necessary information to include in the VAR system that examines commodity price interlinkages over time. Further, across all three periods considered, Economic Output and Money & Credit indicators contribute positively to the factor composition. In particular, the positive contribution from the M2 measure of money supply for the US during the crisis period and that from the Fed Funds Rate post crisis are pronounced; hence, the estimated factors strongly reflect the effect of the Quantitative Easing policy adopted by the US central bank. The contribution of the other categories are mixed, with that from bond returns being noteworthy due to their role as a proxy for long-term interest rates, which impact both the cost of investment in increasing production capacity and on holding inventories, as well as on the composition of asset portfolios across a range of investment possibilities (stocks, bonds, commodities, etc.).
Next, using these estimated factors, we fit a sparse $\mathrm{VAR}(2)$ model to the augmented $(\widehat{F}^\top_t,X^\top_t)^\top$ system. The estimated transition matrices are depicted in Figures (ref) to (ref) as networks\footnote{In all three figures, the left panel corresponds to $\widehat{A}^{(1)}$ and the right panel corresponds to $\widehat{A}^{(2)}$. Node sizes are proportional to node weighted degrees. Positive edges are in red and negative edges are in blue. Edges with higher saturation have larger magnitudes.}. It is apparent that the factors play an important role, both as emitters and receivers. The effects from the first lag are generally stronger to that from the second one. In particular, focusing on the first lag, the dominant nodes in the system have shifted over time from (OIL, SOYBEANS, ZINC) pre crisis to (SUGAR, WHEAT, COPPER) during the crisis, then to (OIL, SOYBEANS, RICE) post crisis. Based on node weighted degree, the role of OIL is dominant in both pre- and post-crisis periods, but is much weaker during the crisis.
Another key feature of the interlinkage networks is their increased connectivity during the crisis period, vis-a-vis the pre- and post-crisis periods. The same empirical finding has been noted for stock returns lin2017regularized. Before the global financial crisis of 2008, commodity prices were fast rising primarily due to increased demand from China. Specifically, as Chinese industrial production quadrupled between 2001 and 2011, its consumption of industrial metals (Copper, Zinc, Aluminum, Lead) increased by 330%, while its oil consumption by 98%. This strong demand shock led to a sharp rise in these commodity prices, particularly accentuated beginning in 2006 (the onset of the crisis period considered in our analysis), briefly disrupted with a quick plunge of commodity prices in 2008 and their subsequent recovery in the ensuing period until late 2010, when demand from China subsided, which coupled with weak demand from the EU, Japan and the US in the aftermath of the crisis created an oversupply that put downward pressure on prices. These events induce strong inter-temporal and cross-temporal correlations amongst commodity prices, and hence are reflected in their estimated interlinkage network.
This paper considered the estimation of FAVAR model under the high-dimensional scaling. It introduced an identifiability constraint (IR+Compactness) that is suitable for high-dimensional settings, and when such a constraint is incorporated in the optimization problem based upon the calibration equation, the global optimizer corresponds to model parameter estimates with bounded statistical errors. This development also allows for accurate estimation of the transition matrices of the VAR system, despite the plug-in factor block contains error due to the fact that it is an estimated quantity. Extensive numerical work illustrates the overall good performance of the proposed empirical implementation procedure, but also illustrates that the imposed constraint is not particularly stringent, especially in settings where the coefficient matrix $\Gamma$ of the observed predictor variables in the calibration equation exhibits sufficient level of sparsity.
The key advantage of the FAVAR model is that it can leverage information from a large number of variables, while modeling the cross-temporal dependencies of a smaller number of them that are of primary interest to the analyst.
Recall that the nature of the FAVAR model results in estimating the transition matrix of a VAR system with one block of the observations (factors) being an estimated quantity, rather than conducting the estimation based on observed samples. Similar in flavor problems have been examined in the high-dimensional iid setting loh2012high, as well as low dimensional time series settings; for example, chanda1996asymptotic examine parameter estimation of a univariate autoregressive process with error-in-variables and in more recent work komunjer2014measurement investigate parameter identification of VAR-X and dynamic panel VAR models subject to measurement errors.
\vskip 0.2in