pid.dc: Independence-based identification of panel SVAR models using distance covariance (DC) statistic

Description

Given an estimated panel of VAR models, this function applies independence-based identification for the structural impact matrix $B_i$ of the corresponding SVAR model $$y_{it} = c_{it} + A_{i1} y_{i,t-1} + ... + A_{i,p_i} y_{i,t-p_i} + u_{it}$$ $$ = c_{it} + A_{i1} y_{i,t-1} + ... + A_{i,p_i} y_{i,t-p_i} + B_i \epsilon_{it}.$$ Matrix $B_i$ corresponds to the unique decomposition of the least squares covariance matrix $\Sigma_{u,i} = B_i B_i'$ if the vector of structural shocks $\epsilon_{it}$ contains at most one Gaussian shock (Comon, 1994). A nonparametric dependence measure, the distance covariance (Szekely et al., 2007), determines least dependent structural shocks. The algorithm described in Matteson and Tsay (2013) is applied to calculate the matrix $B_i$.

Usage

pid.dc(
  x,
  combine = c("group", "pool", "indiv"),
  n.factors = NULL,
  n.iterations = 100,
  PIT = FALSE
)

Value

List of class 'pid' with elements:

A: Matrix. The lined-up coefficient matrices $A_j, j=1,\ldots,p$ for the lagged variables in the panel VAR.
B: Matrix. Mean group of the estimated structural impact matrices $B_i$, i.e. the unique decomposition of the covariance matrices of reduced-form errors.
L.varx: List of 'varx' objects for the individual estimation results to which the structural impact matrices $B_i$ have been added.
eps0: Matrix. The combined whitened residuals $\epsilon_{0}$ to which the ICA procedure has been applied subsequently. These are still the unrotated baseline shocks! NULL if 'indiv' identifications have been used.
ICA: List of objects resulting from the underlying ICA procedure. NULL if 'indiv' identifications have been used.
rotation_angles: Numeric vector. The rotation angles suggested by the combined identification procedure. NULL if 'indiv' identifications have been used.
args_pid: List of characters and integers indicating the identification methods and specifications that have been used.
args_pvarx: List of characters and integers indicating the estimator and specifications that have been used.

Arguments

x: An object of class 'pvarx' or a list of VAR objects that will be coerced to 'varx'. Estimated panel of VAR objects.
combine: Character. The combination of the individual reduced-form residuals via 'group' for the group ICA by Calhoun et al. (2001) using common structural shocks, via 'pool' for the pooled shocks by Herwartz and Wang (2024) using a common rotation matrix, or via 'indiv' for individual-specific $B_i \forall i$ using strictly separated identification runs.
n.factors: Integer. Number of common structural shocks across all individuals if the group ICA is selected.
n.iterations: Integer. The maximum number of iterations in the 'steadyICA' algorithm. The default in 'steadyICA' is 100.
PIT: Logical. If PIT='TRUE', the distribution and density of the independent components are estimated using Gaussian kernel density estimates.

Notes on the Reproduction in "Examples"

The reproduction of Herwartz and Wang (HW, 2024:630) serves as an exemplary application and unit-test of the implementation by pvars. While vars' VAR employs equation-wise lm with the QR-decomposition of the regressor matrix $X$, HW2024 and accordingly the reproduction by pvarx.VAR both calculate $X'(XX')^{-1}$ for the multivariate least-squares estimation of $A_i$. Moreover, both use steadyICA for identification such that the reproduction result for the pooled rotation matrix Q is close to HW2024, the mean absolute difference between both 4x4 matrices is less than 0.0032. Note that the single EA-Model is estimated and identified the same way, which can be extracted as a separate 'varx' object from the trivial panel object by $L.varx[[1]] and even bootstrapped by sboot.mb.

Some differences remain such that the example does not exactly reproduce the results in HW2024. To account for the $n$ exogenous and deterministic regressors in slot $D, pvarx.VAR calculates $\Sigma_{u,i}$ with the degrees of freedom $T-Kp_i-n$ instead of HW2024's $T-Kp_i-1$. Moreover, the confidence bands for the IRF are based on pvars' panel moving-block- instead of HW2024's wild bootstrap. The responses of real GDP and of inflation are not scaled by 0.01, unlike in HW2024. Note that both bootstrap procedures keep D fixed over their iterations.

References

Calhoun, V. D., Adali, T., Pearlson, G. D., and Pekar, J. J. (2001): "A Method for Making Group Inferences from Functional MRI Data using Independent Component Analysis", Human Brain Mapping, 16, pp. 673-690.

Comon, P. (1994): "Independent Component Analysis: A new Concept?", Signal Processing, 36, pp. 287-314.

Herwartz, H., and Wang, S. (2024): "Statistical Identification in Panel Structural Vector Autoregressive Models based on Independence Criteria", Journal of Applied Econometrics, 39 (4), pp. 620-639.

Matteson, D. S., and Tsay, R. S. (2017): "Independent Component Analysis via Distance Covariance", Journal of the American Statistical Association, 112, pp. 623-637.

Risk, B., Matteson, D. S., Ruppert, D., Eloyan, A., and Caffo, B. S. (2014): "An Evaluation of Independent Component Analyses with an Application to Resting-State fMRI", Biometrics, 70, pp. 224-236.

Szekely, G. J., Rizzo, M. L., and Bakirov, N. K. (2007): "Measuring and Testing Dependence by Correlation of Distances", Annals of Statistics, 35, pp. 2769-2794.

Examples

Run this code

### replicate Herwartz,Wang 2024:630, Ch.4 ###
# \donttest{
# select minimal or full example #
is_min = TRUE
n.boot = ifelse(is_min, 5, 1000)
idx_i  = ifelse(is_min, 1, 1:14)

# load and prepare data #
data("EURO")
names_i = names(EURO[idx_i+1])  # country names (#1 is EA-wide aggregated data)
names_s = paste0("epsilon[ ", c(1:2, "m", "f"), " ]")  # shock names
idx_k   = 1:4   # endogenous variables in individual data matrices
idx_t   = 1:(nrow(EURO[[1]])-1)  # periods from 2001Q1 to 2019Q4 
trend2  = idx_t^2

# panel SVARX model, Ch.4.1 #
L.data = lapply(EURO[idx_i+1], FUN=function(x) x[idx_t, idx_k])
L.exog = lapply(EURO[idx_i+1], FUN=function(x) cbind(trend2, x[idx_t, 5:10]))
R.lags = c(1,2,1,2,2,2,2,2,1,2,2,2,2,1)[idx_i]; names(R.lags) = names_i
R.pvar = pvarx.VAR(L.data, lags=R.lags, type="both", D=L.exog)
R.pid  = pid.dc(R.pvar, combine="pool")
print(R.pid)  # suggests e3 and e4 to be MP and financial shocks, respectively.
colnames(R.pid$B) = names_s  # accordant labeling

# EA-wide SVARX model, Ch.4.2 #
R.data = EURO[[1]][idx_t, idx_k]
R.exog = cbind(trend2, EURO[[1]][idx_t, 5:6])
R.varx = pvarx.VAR(list(EA=R.data), lags=2, type="both", D=list(EA=R.exog))
R.id   = pid.dc(R.varx, combine="indiv")$L.varx$EA
colnames(R.id$B) = names_s  # labeling

# comparison of IRF without confidence bands, Ch.4.3.1 #
data("EU_w")  # GDP weights with the same ordering names_i as L.varx in R.pid
R.norm = function(B) B / matrix(diag(B), nrow(B), ncol(B), byrow=TRUE) * 25
R.irf  = as.pplot(
  EA=plot(irf(R.id,  normf=R.norm), selection=list(idx_k, 3:4)),
  MG=plot(irf(R.pid, normf=R.norm, w=EU_w), selection=list(idx_k, 3:4)),
  color_g=c("#3B4992FF", "#008B45FF"), shape_g=16:17, n.rows=length(idx_k))
plot(R.irf)

# comparison of IRF with confidence bands, Ch.4.3.1 #
R.boot_EA = sboot.mb(R.id, b.length=8, n.boot=n.boot, n.cores=2, normf=R.norm)
R.boot_MG = sboot.pmb(R.pid, b.dim=c(8, FALSE), n.boot=n.boot, n.cores=2, 
                      normf=R.norm, w=EU_w)
R.irf = as.pplot(
  EA=plot(R.boot_EA, lowerq=0.16, upperq=0.84, selection=list(idx_k, 3:4)),
  MG=plot(R.boot_MG, lowerq=0.16, upperq=0.84, selection=list(idx_k, 3:4)),
  color_g=c("#3B4992FF", "#008B45FF"), shape_g=16:17, n.rows=length(idx_k))
plot(R.irf)
# }

Run the code above in your browser using DataLab