Package 'pvars'

Coerce into a 'pplot' object

Description

Coerce into a 'pplot' object. This function allows (1) to overlay and colorize multiple plots of IRF and confidence bands in a single 'ggplot' object and (2) to overwrite shock and variable names in verbatim LaTeX math mode suitable for the export via tikzDevice.

Usage

as.pplot(
  ...,
  names_k = NULL,
  names_s = NULL,
  names_g = NULL,
  color_g = NULL,
  shape_g = NULL,
  n.rows = NULL,
  scales = "free_y",
  Latex = FALSE
)

Arguments

...	A single ggplot or list(s) of ggplots to be transformed.
names_k	Vector. Names of the variables $k=1,\ldots,K$. If NULL (the default), the names of the first plot are reused. For LaTeX exports, use e.g. paste0("y_{ ", 1:dim_K, " }").
names_s	Vector. Names of the shocks $s=1,\ldots,S$. If NULL (the default), the names of the first plot are reused. For LaTeX exports, use e.g. paste0("\\epsilon_{ ", 1:dim_S, " }").
names_g	Vector. Names of the layered plots $g=1,\ldots,G$. If NULL (the default), the names of the plots given in ... are reused.
color_g	Vector. Colors of the layered plots $g=1,\ldots,G$.
shape_g	Vector. Shapes of the layered plots $g=1,\ldots,G$, see ggplot2's geom_point for shapes. If NULL (the default), no points will be set on the IRF-lines.
n.rows	Integer. Number of rows in facet_wrap. If NULL (the default), the dimensions of the facet plots given in ... are reused.
scales	Character. Should scales be fixed ("fixed"), free ("free"), or free in one dimension ("free_x", "free_y", the default)? See facet_wrap.
Latex	Logical. If TRUE, the arrows of the facet labels are written in LaTeX math mode.

Details

as.pplot is used as an intermediary in the 'pplot' functions to achieve compatibility with different classes. Equivalently to as.pvarx for panels of $N$ VAR objects, as.pplot summarizes panels of $G$ plot objects that have been returned from the 'plot' method for class 'svarirf' or 'sboot'. If the user wishes to extend the compatibility of the 'pplot' functions with further classes, she may simply specify accordant as.pplot-methods instead of altering the original 'pplot' functions.

Value

A list of class 'pplot'. Objects of this class contain the elements:

F.plot	'ggplot' object for the merged plot.
L.plot	List of 'ggplot' objects containing all $G$ plots.
args_pplot	List of characters and integers indicating the specifications used for creating F.plot.

See Also

PP, irf.pvarx, pid.dc, and id.iv for further examples of edited plots, in particular for subset and reordered facet plots with reshaped facet dimensions.

Examples

### gallery of merged IRF plots ###

library("ggplot2")
data("PCAP")
names_k = c("g", "k", "l", "y")  # variable names
names_i = levels(PCAP$id_i)      # country names 
L.data  = sapply(names_i, FUN=function(i) 
  ts(PCAP[PCAP$id_i==i, names_k], start=1960, end=2019, frequency=1), 
  simplify=FALSE)
L.vars = lapply(L.data, FUN=function(x) vars::VAR(x, p=2, type="both"))
L.chol = lapply(L.vars, FUN=function(x) svars::id.chol(x))

# overlay all IRF to get an overview on the stability #
L.irf = lapply(L.chol, FUN=function(x) plot(irf(x, n.ahead=30)))
summary(as.pvarx(L.vars))
as.pplot(L.irf)

# overlay IRF of selected countries and quantiles of all countries #
F.mg  = plot(sboot.mg(L.chol, n.ahead=30), lowerq=0.05, upperq=0.95)
R.irf = as.pplot(MG=F.mg, L.irf[c("DEU", "FRA", "ITA", "JPN")])
plot(R.irf)  # emphasize MG-IRF in next step
R.irf = as.pplot(R.irf, color_g="black", shape_g=c(20, rep(NA, 4)))
R.irf$F.plot + guides(fill="none") + labs(color="Country", shape="Country")

# compare two mean-groups and their quantiles #
idx_nord = c(5, 6, 10, 17, 20)  # Nordic countries
R.irf = as.pplot(color_g=c("black", "blue"), 
  Other  = plot(sboot.mg(L.chol[-idx_nord])), 
  Nordic = plot(sboot.mg(L.chol[ idx_nord])))
plot(R.irf)
 
# compare different shock ordering for MG-IRF #
R.pid1 = pid.chol(L.vars)
R.pid2 = pid.chol(L.vars, order_k=4:1)
R.pid3 = pid.chol(L.vars, order_k=c(1,4,2,3))

R.pal = RColorBrewer::brewer.pal(n=8, name="Spectral")[c(8, 1, 4)]
R.irf = as.pplot(color_g=R.pal, shape_g=c(2, 3, 20), 
  GKLY = plot(irf(R.pid1, n.ahead=25)), 
  YLKG = plot(irf(R.pid2, n.ahead=25)), 
  GYKL = plot(irf(R.pid3, n.ahead=25)))
R.mg = as.pplot(color_g=R.pal, shape_g=c(2, 3, 20), 
  GKLY = plot(sboot.mg(R.pid1, n.ahead=25), lowerq=0.05, upperq=0.95), 
  YLKG = plot(sboot.mg(R.pid2, n.ahead=25), lowerq=0.05, upperq=0.95), 
  GYKL = plot(sboot.mg(R.pid3, n.ahead=25), lowerq=0.05, upperq=0.95))
         
# colorize and export a single sub-plot to Latex #
library("tikzDevice")
textwidth = 15.5/2.54  # LaTeX textwidth from "cm" into "inch"
file_fig  = file.path(tempdir(), "Fig_irf.tex")

R.irf = as.pplot(
  DEU = plot(irf(L.chol[["DEU"]], n.ahead=50), selection=list(4, 1)), 
  FRA = plot(irf(L.chol[["FRA"]], n.ahead=50), selection=list(4, 1)),
  color_g = c("black", "blue"),
  names_g = c("Germany", "France"),
  names_k = "y", 
  names_s = "\\epsilon_{ g }", 
  Latex   = TRUE)

tikz(file=file_fig, width=1.2*textwidth, height=0.8*textwidth)
  R.irf$F.plot + labs(color="Country") + theme_minimal()
dev.off()

---

Coerce into a 'pvarx' object

Description

Coerce into a 'pvarx' object. On top of the parent class 'pvarx', the child class 'pid' is imposed if the input object to be transformed contains a panel SVAR model.

Usage

as.pvarx(x, w = NULL, ...)

Arguments

x	A panel VAR object to be transformed.
w	Numeric, logical, or character vector. $N$ numeric elements weighting the individual coefficients, or names or $N$ logical elements selecting a subset from the individuals $i = 1, \ldots, N$ for the MG estimation. If NULL (the default), all $N$ individuals are included without weights.
...	Additional arguments to be passed to or from methods.

Details

as.pvarx is used as an intermediary in the pvars functions to achieve compatibility with different classes of panel VAR objects. If the user wishes to extend this compatibility with further classes, she may simply specify accordant as.pvarx-methods instead of altering the original pvars function.

Value

A list of class 'pvarx'. Objects of this class contain the elements:

A	Matrix. The lined-up coefficient matrices $A_j, j=1,\ldots,p$ for the lagged variables in the panel VAR.
B	Matrix. The $(K \times S)$ structural impact matrix of the panel SVAR model or an identity matrix $I_K$ as a placeholder for the unidentified VAR model.
beta	Matrix. The $((K+n_{d1}) \times r)$ cointegrating matrix of the VAR model if transformed from a rank-restricted VECM.
L.varx	List of varx objects for the individual estimation results.
args_pvarx	List of characters and integers indicating the estimator and specifications that have been used.
args_pid	List of characters and integers indicating the identification methods and specifications that have been used. This element is specific to the child-class 'pid' for panel SVAR models, that inherit from parent-class 'pvarx' for any panel VAR model.

Examples

data("PCAP")
names_k = c("g", "k", "l", "y")  # variable names
names_i = levels(PCAP$id_i)      # country names 
L.data  = sapply(names_i, FUN=function(i) 
  ts(PCAP[PCAP$id_i==i, names_k], start=1960, end=2019, frequency=1), 
  simplify=FALSE)
  
L.vars = lapply(L.data, FUN=function(x) vars::VAR(x, p=2, type="both"))
as.pvarx(L.vars)

---

Deterministic regressors in pvars

Description

Deterministic regressors can be specified via the arguments of the conventional 'type', customized 'D', and period-specific 't_D'. While 'type' is a single character and 'D' a data matrix of dimension $(n_{\bullet} \times (p+T))$, the specifications for $\tau$ in the list 't_D' are more complex and therefore preventively checked by as.t_D.

Usage

as.t_D(x, ...)

Arguments

x	A list of vectors for $\tau$ to be checked. Since 'x' is just checked, Section "Value" explains function-input and -output likewise.
...	Additional arguments to be passed to or from methods.

Value

A list of class 't_D' specifying $\tau$. Objects of this class can exclusively contain the elements:

t_break	Vector of integers. The activating periods for trend breaks $d = [\ldots, 0, 0, 1, 2, 3, \ldots$].
t_shift	Vector of integers. The activating periods for shifts in the constant $d = [\ldots, 0, 0, 1, 1, 1, \ldots$].
t_impulse	Vector of integers. The activating periods for single impulses $d = [\ldots, 0, 0, 1, 0, 0, \ldots$].
t_blip	Vector of integers. The activating period for blips $d = [\ldots, 0, 0, 1, -1, 0, \ldots$].
n.season	Integer. The number of seasons.

Reference Time Interval

The complete time series (i.e. including the presample) constitutes the reference time interval. Accordingly, 'D' contains $p+T$ observations, and 't_D' contains the positions of activating periods $\tau$ in $1,\ldots,(p+T)$. In a balanced panel $p_i+T_i = T^*$, the same date implies the same $\tau$ in $1,\ldots,T^*$, as shown in the example for pcoint.CAIN. However, in an unbalanced panel, the same date can imply different $\tau$ across $i$ in accordance with the individual time interval $1,\ldots,(p_i+T_i)$. Note that across the time series in 'L.data', it is the last observation in each data matrix that refers to the same date.

Conventional Type

An overview is given here and a detailed explanation in the package vignette.

• type (VAR) is specified in VAR models just as in vars' VAR, namely by a 'const', a linear 'trend', 'both', or 'none' of those.

• type_SL is used in the 'additive' SL procedure for testing the cointegration rank only, which removes the mean ('SL_mean') or mean and linear trend ('SL_trend') by GLS-detrending.

• type (VECM) is used in the 'innovative' Johansen procedure for testing the cointegration rank and estimating the VECM. In accordance with Juselius (2007, Ch.6.3), the available model specifications are: 'Case1' for none, 'Case2' for a constant in the cointegration relation, 'Case3' for an unrestricted constant, 'Case4' for a linear trend in the cointegration relation and an unrestricted constant, or 'Case5' for an unrestricted constant and linear trend.

References

Juselius, K. (2007): The Cointegrated VAR Model: Methodology and Applications, Advanced Texts in Econometrics, Oxford University Press, USA, 2nd ed.

Examples

t_D = list(t_impulse=c(10, 20, 35), t_shift=10)
as.t_D(t_D)

---

Coerce into a 'varx' object

Description

Coerce into a 'varx' object. On top of the parent class 'varx', the child class 'id' is imposed if the input object to be transformed contains an SVAR model.

Usage

as.varx(x, ...)

Arguments

x	A VAR object to be transformed.
...	Additional arguments to be passed to or from methods.

Details

as.varx is used as an intermediary in the pvars functions to achieve compatibility with different classes of VAR objects. If the user wishes to extend this compatibility with further classes, she may simply specify accordant as.varx-methods instead of altering the original pvars function. Classes already covered by pvars are those of the vars ecosystem, in particular the classes

• 'varest' for reduced-form VAR estimates from VAR,

• 'vec2var' for reduced-form VECM estimates from vec2var,

• 'svarest' for structural VAR estimates from BQ,

• 'svecest' for structural VECM estimates from SVEC, and

• 'svars' for structural VAR estimates from svars' id.chol, id.cvm, or id.dc.

By transformation to 'varx', these VAR estimates can thus be subjected to pvars' bootstrap procedure sboot.mb and S3 methods such as summary and toLatex.

Value

A list of class 'varx'. Objects of this class contain the elements:

A	Matrix. The lined-up VAR coefficient matrices $A_j, j=1,\ldots,p$ for the lagged variables.
B	Matrix. The $(K \times S)$ structural impact matrix of the SVAR model or an identity matrix $I_K$ as a placeholder for the unidentified VAR model.
beta	Matrix. The $((K+n_{d1}) \times r)$ cointegrating matrix of the VAR model if transformed from a rank-restricted VECM.
SIGMA	Matrix. The $(K \times K)$ residual covariance matrix estimated by least-squares.

The following integers indicate the size of dimensions:

dim_K	Integer. The number of endogenous variables $K$ in the full-system.
dim_S	Integer. The number of identified shocks $S$ in the SVAR model.
dim_T	Integer. The number of time periods $T$ without presample.
dim_p	Integer. The lag-order $p$ of the VAR model in levels.
dim_r	Integer. The cointegration rank $r$ of the VAR model if transformed from a rank-restricted VECM.

Some further elements required for the bootstrap functions are:

y	Matrix. The $(K \times (p+T))$ endogenous variables.
D, D1, D2	Matrices. The $(n_{\bullet} \times (p+T))$ deterministic variables, fixed over bootstrap iterations, (un)restricted to the cointegration relations of the VAR model if transformed from a rank-restricted VECM.
resid	Matrix. The $(K \times T)$ residual matrix.
args_varx	List of characters and integers indicating the estimator and specifications that have been used.
args_id	List of characters and integers indicating the identification methods and specifications that have been used. This element is specific to the child-class 'id' for SVAR models, that inherit from parent-class 'varx' for any VAR model.

Examples

data("PCIT")
names_k = c("APITR", "ACITR", "PITB", "CITB", "GOV", "RGDP", "DEBT")

# estimate reduced-form VAR and coerce into 'varx' object #
R.vars = vars::VAR(PCIT[ , names_k], p=4, type="const")
as.varx(R.vars)

# identify structural VAR and coerce into 'id' object #
R.svar = svars::id.chol(R.vars, order_k=names_k)
as.varx(R.svar)

---

Test procedures for the cointegration rank

Description

Performs test procedures for the rank of cointegration in a single VAR model. The $p$-values are approximated by gamma distributions, whose moments are automatically adjusted to potential period-specific deterministic regressors and weakly exogenous regressors in the partial VECM.

Usage

coint.JO(
  y,
  dim_p,
  x = NULL,
  dim_q = dim_p,
  type = c("Case1", "Case2", "Case3", "Case4", "Case5"),
  t_D1 = NULL,
  t_D2 = NULL
)

coint.SL(y, dim_p, type_SL = c("SL_mean", "SL_trend"), t_D = NULL)

Arguments

y	Matrix. A $(K \times (p+T))$ data matrix of the $K$ endogenous time series variables.
dim_p	Integer. Lag-order $p$ for the endogenous variables y.
x	Matrix. A $(L \times (q+T))$ data matrix of the $L$ weakly exogenous time series variables.
dim_q	Integer. Lag-order $q$ for the weakly exogenous variables x. The literature uses dim_p (the default).
type	Character. The conventional case of the deterministic term in the Johansen procedure.
t_D1	List of vectors. The activating break periods $\tau$ for the period-specific deterministic regressors in $d_{1,t}$, which are restricted to the cointegration relations. The accompanying lagged regressors are automatically included in $d_{2,t}$. The $p$-values are calculated for up to two breaks resp. three sub-samples.
t_D2	List of vectors. The activating break periods $\tau$ for the period-specific deterministic regressors in $d_{2,t}$, which are unrestricted.
type_SL	Character. The conventional case of the deterministic term in the Saikkonen-Luetkepohl (SL) procedure.
t_D	List of vectors. The activation periods $\tau$ for the period-specific deterministic regressors in $d_{t}$ of the SL-procedure. The accompanying lagged regressors are automatically included in $d_{t}$. The $p$-values are calculated for up to two breaks resp. three sub-samples.

Value

A list of class 'coint', which contains elements of length $K$ for each $r_{H0}=0,\ldots,K-1$:

r_H0	Rank under each null hypothesis.
stats_TR	Trace (TR) test statistics.
stats_ME	Maximum eigenvalue (ME) test statistics.
pvals_TR	$p$-values of the TR test.
pvals_ME	$p$-values of the ME test. NA if moments of the gamma distribution are not available for the chosen data generating process.
lambda	Eigenvalues, the squared canonical correlation coeffcients (saved only for the Johansen procedure).
args_coint	List of characters and integers indicating the cointegration test and specifications that have been used.

Functions

• coint.JO(): Johansen procedure.

• coint.SL(): (Trenkler)-Saikkonen-Luetkepohl procedure.

References

Johansen, S. (1988): "Statistical Analysis of Cointegration Vectors", Journal of Economic Dynamics and Control, 12, pp. 231-254.

Doornik, J. (1998): "Approximations to the Asymptotic Distributions of Cointegration Tests", Journal of Economic Surveys, 12, pp. 573-93.

Johansen, S., Mosconi, R., and Nielsen, B. (2000): "Cointegration Analysis in the Presence of Structural Breaks in the Deterministic Trend", Econometrics Journal, 3, pp. 216-249.

Kurita, T., Nielsen, B. (2019): "Partial Cointegrated Vector Autoregressive Models with Structural Breaks in Deterministic Terms", Econometrics, 7, pp. 1-35.

Saikkonen, P., and Luetkepohl, H. (2000): "Trend Adjustment Prior to Testing for the Cointegrating Rank of a Vector Autoregressive Process", Journal of Time Series Analysis, 21, pp. 435-456.

Trenkler, C. (2008): "Determining $p$-Values for Systems Cointegration Tests with a Prior Adjustment for Deterministic Terms", Computational Statistics, 23, pp. 19-39.

Trenkler, C., Saikkonen, P., and Luetkepohl, H. (2008): "Testing for the Cointegrating Rank of a VAR Process with Level Shift and Trend Break", Journal of Time Series Analysis, 29, pp. 331-358.

Examples

### reproduce basic example in "urca" ###
library("urca")
data(denmark)
sjd = denmark[ , c("LRM", "LRY", "IBO", "IDE")]

# rank test and estimation of the full VECM as in "urca" #
R.JOrank = coint.JO(y=sjd,      dim_p=2, type="Case2", t_D2=list(n.season=4))
R.JOvecm = VECM(y=sjd, dim_r=1, dim_p=2, type="Case2", t_D2=list(n.season=4))

# ... and of the partial VECM, i.e. after imposing weak exogeneity #
R.KNrank = coint.JO(y=sjd[ , c("LRM"), drop=FALSE],   dim_p=2,
                    x=sjd[ , c("LRY", "IBO", "IDE")], dim_q=2, 
                    type="Case2", t_D1=list(t_shift=36), t_D2=list(n.season=4))
R.KNvecm = VECM(y=sjd[ , c("LRM"), drop=FALSE],   dim_p=2,
                x=sjd[ , c("LRY", "IBO", "IDE")], dim_q=2, dim_r=1, 
                type="Case2", t_D1=list(t_shift=36), t_D2=list(n.season=4))

### reproduce Oersal,Arsova 2016:22, Tab.7.5 "France" ###
data("ERPT")
names_k = c("lpm5", "lfp5", "llcusd")  # variable names for "Chemicals and related products"
names_i = levels(ERPT$id_i)[c(1,6,2,5,4,3,7)]  # ordered country names
L.data  = sapply(names_i, FUN=function(i) 
   ts(ERPT[ERPT$id_i==i, names_k], start=c(1995, 1), frequency=12), 
   simplify=FALSE)
R.TSLrank = coint.SL(y=L.data$France, dim_p=3, type_SL="SL_trend", t_D=list(t_break=89))

---

Data set on the Exchange Rate Pass-Through

Description

The data set ERPT consists of monthly observations for the logarithm of import prices $lm^*$, foreign prices $lf^*$, and the exchange rate against the US dollar $llcusd$. It covers the period January 1995 to March 2005 $(T=123)$ for $N=7$ countries. The asterisk denotes the industry of the variables and can take values from 0 to 8:

• 0: Food and live animals chiefly for food

• 1: Beverages and tobacco

• 2: Crude materials (inedible, except fuels)

• 3: Mineral fuels, lubricants and related materials

• 4: Animal and vegetable oils, fats and waxes

• 5: Chemicals and related products

• 6: Manufactured goods classified chiefly by materials

• 7: Machines, transport equipment

• 8: Manufactured goods

Usage

data("ERPT")

Format

A long-format data panel of class 'data.frame', where the columns id_i and id_t indicate the country and month respectively.

Source

The prepared Eurostat data set is directly obtainable from the ZBW Journal Data Archive: doi:10.15456/jae.2022321.0717881037. This is open data under the CC BY 4.0 license.

References

Banerjee, A., and Carrion-i-Silvestre, J. L. (2015): "Cointegration in Panel Data with Structural Breaks and Cross-Section Dependence", Journal of Applied Econometrics, 30 (1), pp. 1-23.

See Also

Other data sets: EURO, EU_w, ICAP, MDEM, MERM, PCAP, PCIT

---

Weights for the Euro Monetary Policy Transmission

Description

The data set EU_w is a vector of 14 elements. These are weights for $N=14$ member countries of the Euro area, constructed as the average share of their respective real GDP over the sample period in Herwartz, Wang (2024).

Usage

data("EU_w")

Format

A numeric vector containing 14 named elements.

Source

The prepared Eurostat data set is directly obtainable from the ZBW Journal Data Archive: doi:10.15456/jae.2024044.1425287131. This is open data under the CC BY 4.0 license in accordance with the deposit license of the ZBW Journal Data Archive.

References

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.

See Also

Other data sets: ERPT, EURO, ICAP, MDEM, MERM, PCAP, PCIT

---

Data set on the Euro Monetary Policy Transmission

Description

The data set EURO is a list of 15 'data.frame' objects, each consisting of quarterly observations for

• the first-difference of log real GDP on national $dl\_GDP$ or aggregated EA-level $EA\_dl\_GDP$,

• the annualized inflation of the (log) GDP deflator on national $dl\_deflator$ or aggregated EA-level $EA\_pi$,

• the EA-wide short-term interest rate $IR$,

• the EA-wide option-adjusted bond spreads $BBB$,

• the first-difference of log real GDP in the remaining countries $dl\_GDP\_EA$,

• the weighted inflation in the remaining countries $dl\_deflator\_EA$,

• the inflation of a world commodity price index $WCP$,

• the US effective federal funds rate $US\_FFR$,

• the trade volume in percentage of GDP $trade$, and

• the government spending in percentage of GDP $ge$.

The data covers the period Q1 2001 to Q1 2020 $(T=77)$ for the aggregate of the Euro area (EA, first element in list) and $N=14$ of its member countries (subsequent 14 elements in list).

Usage

data("EURO")

Format

A list-format data panel of class 'list' containing 15 'data.frame' objects with named time series.

Source

The prepared Eurostat data set is directly obtainable from the ZBW Journal Data Archive: doi:10.15456/jae.2024044.1425287131. This is open data under the CC BY 4.0 license in accordance with the deposit license of the ZBW Journal Data Archive.

References

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.

See Also

Other data sets: ERPT, EU_w, ICAP, MDEM, MERM, PCAP, PCIT

---

Forecast Error Variance Decomposition

Description

Calculates the forecast error variance decomposition. Respects SVAR models of cases $S \neq K$, i.e. partially identified or excess shocks, too.

Usage

## S3 method for class 'id'
fevd(x, n.ahead = 10, ...)

Arguments

x	SVAR object of class 'id' or any other that will be coerced to ('id', 'varx').
n.ahead	Integer specifying the steps ahead, i.e. the horizon of the FEVD.
...	Currently not used.

Value

A list of class 'svarfevd' holding the forecast error variance decomposition of each variables as a 'data.frame'.

References

Luetkepohl, H. (2005): New Introduction to Multiple Time Series Analysis, Springer, 2nd ed.

Jentsch, Lunsford (2022): "Asymptotically Valid Bootstrap Inference for Proxy SVARs", Journal of Business and Economic Statistics, 40, pp. 1876-1891.

Examples

data("PCIT")
names_k = c("APITR", "ACITR", "PITB", "CITB", "GOV", "RGDP", "DEBT")
names_l = c("m_PI", "m_CI")  # proxy names
names_s = paste0("epsilon[ ", c("PI", "CI"), " ]")  # shock names
dim_p   = 4  # lag-order

# estimate and identify proxy SVAR #
R.vars = vars::VAR(PCIT[ , names_k], p=dim_p, type="const")
R.idBL = id.iv(R.vars, iv=PCIT[-(1:dim_p), names_l], S2="MR", cov_u="OMEGA")
colnames(R.idBL$B) = names_s  # labeling

# calculate and plot FEVD under partial identification #
plot(fevd(R.idBL, n.ahead=20))

---

Data set on Infrastructure Capital Stocks

Description

The data set ICAP consists of annual observations for

• the real GDP $Y$ in US dollars,

• the physical capital stocks $K$ in US dollars,

• the physical capital stocks $KBC30$ via "backcasting" in US dollars,

• the number of workers $LWDI$ providing the total labor force, and

• the average years of secondary education $secondary$ of the population.

It further reports physical measures of infrastructure given by

• the electricity generation capacity $EGC$ in megawatts,

• the number of main phone lines $mlines$,

• the kilometers of total roads $troads$,

• the number of cell phones lines $cells$,

• the kilometers of paved roads $proads$, and

• the kilometers of rails $rails$.

It covers the period 1960 to 2000 $(T=41)$ for $N=97$ countries. The monetary values are given in US-Dollars at 2000 prices, i.e. constant PPP.

Usage

data("ICAP")

Format

A long-format data panel of class 'data.frame', where the columns id_i and id_t indicate the country and year respectively. Column COUNTRY contains the complete country names.

Source

The prepared data set is directly obtainable from the ZBW Journal Data Archive: doi:10.15456/jae.2022321.0717368489. This is open data under the CC BY 4.0 license.

References

Calderon, C., Moral-Benito, E., and Serven, L. (2015): "Is Infrastructure Capital Productive? A Dynamic Heterogeneous Approach", Journal of Applied Econometrics, 30 (2), pp. 177-198.

See Also

Other data sets: ERPT, EURO, EU_w, MDEM, MERM, PCAP, PCIT

---

Identification of SVEC models by imposing long- and short-run restrictions

Description

Identifies an SVEC model by utilizing a scoring algorithm to impose long- and short-run restrictions. See the details of SVEC in vars.

Usage

id.grt(
  x,
  LR = NULL,
  SR = NULL,
  start = NULL,
  max.iter = 100,
  conv.crit = 1e-07,
  maxls = 1
)

Arguments

x	VAR object of class 'varx' estimated under rank-restriction or any other that will be coerced to 'varx'.
LR	Matrix. The restricted long-run impact matrix.
SR	Matrix. The restricted contemporaneous impact matrix.
start	Vector. The starting values for $\gamma$, set by rnorm if NULL (the default).
max.iter	Integer. The maximum number of iterations.
conv.crit	Real number. Convergence value of algorithm.
maxls	Real number. Maximum movement of the parameters between two iterations of the scoring algorithm.

Value

List of class 'id'.

References

Amisano, G. and Giannini, C. (1997): Topics in Structural VAR Econometrics, Springer, 2nd ed.

Breitung, J., Brueggemann R., and Luetkepohl, H. (2004): "Structural Vector Autoregressive Modeling and Impulse Responses", in Applied Time Series Econometrics, ed. by H. Luetkepohl and M. Kraetzig, Cambridge University Press, Cambridge.

Johansen, S. (1996): Likelihood-Based Inference in Cointegrated Vector Autoregressive Models, Advanced Texts in Econometrics, Oxford University Press, USA.

Luetkepohl, H. (2005): New Introduction to Multiple Time Series Analysis, Springer, 2nd ed.

Pfaff, B. (2008): "VAR, SVAR and SVEC Models: Implementation within R Package vars", Journal of Statistical Software, 27, pp. 1-32.

See Also

... the original SVEC by Pfaff (2008) in vars. Note that id.grt is just a graftage, but allows for the additional model specifications in VECM and for the bootstrap procedures in sboot.mb, both provided by the pvars package.

Other identification functions: id.iv()

Examples

### reproduce basic example in "vars" ###
library(vars)
data("Canada")
names_k = c("prod", "e", "U", "rw")  # variable names
names_s = NULL  # optional shock names

# colnames of the restriction matrices are passed as shock names #
SR = matrix(NA, nrow=4, ncol=4, dimnames=list(names_k, names_s))
SR[4, 2] = 0
LR = matrix(NA, nrow=4, ncol=4, dimnames=list(names_k, names_s))
LR[1, 2:4] = 0
LR[2:4, 4] = 0

# estimate and identify SVECM #
R.vecm = VECM(y=Canada[ , names_k], dim_p=3, dim_r=1, type="Case4")
R.grt  = id.grt(R.vecm, LR=LR, SR=SR)

---

Identification of SVAR models by means of proxy variables

Description

Given an estimated VAR model, this function uses proxy variables to partially identify the structural impact matrix $B$ of the corresponding SVAR model

$y_{t} = c_{t} + A_{1} y_{t-1} + ... + A_{p} y_{t-p} + u_{t}$

$= c_{t} + A_{1} y_{t-1} + ... + A_{p} y_{t-p} + B \epsilon_{t}.$

In general, identification procedures determine $B$ up to column ordering, scale, and sign. For a unique solution, id.iv follows the literature on proxy SVAR. The $S$ columns in $B = [B_1 : B_2]$ of the identified shocks $\epsilon_{ts}, s=1,\ldots,S,$ are ordered first, and the variance $\sigma^2_{\epsilon,s} = 1$ is normalized to unity (see e.g. Lunsford 2015:6, Eq. 9). Further, the sign is fixed to a positive correlation between proxy and shock series. A normalization of the impulsed shock that may fix the size of the impact response in the IRF can be imposed subsequently via 'normf' in irf.varx and sboot.mb.

Usage

id.iv(x, iv, S2 = c("MR", "JL", "NQ"), cov_u = "OMEGA", R0 = NULL)

Arguments

x	VAR object of class 'varx' or any other that will be coerced to 'varx'.
iv	Matrix. A $(L \times T)$ data matrix of the $L$ proxy time series $m_t$.
S2	Character. Identification within multiple proxies $m_t$ via 'MR' for lower-triangular $[I_S : -B_{11} B_{12}^{-1} ] B_{1}$ by Mertens, Ravn (2013), via 'JL' for chol$(\Sigma_{mu} \Sigma_{u}^{-1} \Sigma_{um})$ by Jentsch, Lunsford (2021), or via 'NQ' for the nearest orthogonal matrix from svd decomposition by Empting et al. (2025). In case of $S=L=1$ proxy, all three choices provide identical results on $B_1$.
cov_u	Character. Selection of the estimated residual covariance matrix $\hat{\Sigma}_{u}$ to be used in the identification procedure. Either 'OMEGA' (the default) for $\hat{U} \hat{U}'/T_i$ as used in Mertens, Ravn (2013) and Jentsch, Lunsford (2021) or 'SIGMA' for $\hat{U}\hat{U}'/(T-n_{z})$, which corrects for the number of regressors $n_z$. Both character options refer to the name of the respective estimate in the 'varx' object.
R0	Matrix. A $(L \times S)$ selection matrix for 'NQ' that governs the attribution of the $L$ proxies to their specific $S$ structural shock series. If NULL (the default), R0 $= I_S$ will be used such that the $S=L$ columns of $B_1$ are one-by-one estimated from the $S=L$ proxy series 'iv' available.

Value

List of class 'id'.

References

Mertens, K., and Ravn, M. O. (2013): "The Dynamic Effects of Personal and Corporate Income Tax Changes in the United States", American Economic Review, 103, pp. 1212-1247.

Lunsford, K. G. (2015): "Identifying Structural VARs with a Proxy Variable and a Test for a Weak Proxy", Working Paper, No 15-28, Federal Reserve Bank of Cleveland.

Jentsch, C., and Lunsford, K. G. (2019): "The Dynamic Effects of Personal and Corporate Income Tax Changes in the United States: Comment", American Economic Review, 109, pp. 2655-2678.

Jentsch, C., and Lunsford, K. G. (2021): "Asymptotically Valid Bootstrap Inference for Proxy SVARs", Journal of Business and Economic Statistics, 40, pp. 1876-1891.

Empting, L. F. T., Maxand, S., Oeztuerk, S., and Wagner, K. (2025): "Inference in Panel SVARs with Cross-Sectional Dependence of Unknown Form".

See Also

... the individual identification approaches by Lange et al. (2021) in svars.

Other identification functions: id.grt()

Examples

### reproduce Jentsch,Lunsford 2019:2668, Ch.III ###
data("PCIT")
names_k = c("APITR", "ACITR", "PITB", "CITB", "GOV", "RGDP", "DEBT")
names_l = c("m_PI", "m_CI")  # proxy names
names_s = paste0("epsilon[ ", c("PI", "CI"), " ]")  # shock names
dim_p   = 4  # lag-order

# estimate and identify under ordering "BLUE" of Fig.1 and 2 #
R.vars = vars::VAR(PCIT[ , names_k], p=dim_p, type="const")
R.idBL = id.iv(R.vars, iv=PCIT[-(1:dim_p), names_l], S2="MR", cov_u="OMEGA")
colnames(R.idBL$B) = names_s  # labeling

# estimate and identify under ordering "RED" of Fig.1 and 2 #
R.vars = vars::VAR(PCIT[ , names_k[c(2:1, 3:7)]], p=dim_p, type="const")
R.idRD = id.iv(R.vars, iv=PCIT[-(1:dim_p),names_l[2:1]], S2="MR", cov_u="OMEGA")
colnames(R.idRD$B) = names_s[2:1]  # labeling


# select minimal or full example #
is_min = TRUE
n.boot = ifelse(is_min, 5, 10000)

# bootstrap both under 1%-response normalization #
set.seed(2389)
R.norm = function(B) B / matrix(-diag(B), nrow(B), ncol(B), byrow=TRUE)
R.sbBL = sboot.mb(R.idBL, b.length=19, n.boot=n.boot, normf=R.norm)
R.sbRD = sboot.mb(R.idRD, b.length=19, n.boot=n.boot, normf=R.norm)

# plot IRF of Fig.1 and 2 with 68%-confidence levels #
library("ggplot2")
L.idx = list(BLUE1=c(1, 11, 5, 7, 3,  9)+0.1,
             RED1 =c(4, 12, 6, 8, 2, 10)+0.1, 
             RED2 =c(1, 11, 7, 5, 3,  9)+0.1, 
             BLUE2=c(4, 12, 8, 6, 2, 10)+0.1)
# Indexes to subset and reorder sub-plots in plot.sboot(), where 
# the 14 IRF-subplots in the 2D-panel are numbered as a 1D-sequence 
# vectorized by row. '+0.1' makes sub-setting robust against 
# truncation errors from as.integer(). In a given figure, the plots
# RED and BLUE display the same selection of IRF-subplots. 

R.fig1 = as.pplot(
 BLUE=plot(R.sbBL, lowerq=0.16, upperq=0.84, selection=list(1, L.idx[[1]])),
 RED =plot(R.sbRD, lowerq=0.16, upperq=0.84, selection=list(1, L.idx[[2]])),
 names_g=c("APITR first", "ACITR first"), color_g=c("blue", "red"), n.rows=3)

R.fig2 = as.pplot(
 RED =plot(R.sbRD, lowerq=0.16, upperq=0.84, selection=list(1, L.idx[[3]])), 
 BLUE=plot(R.sbBL, lowerq=0.16, upperq=0.84, selection=list(1, L.idx[[4]])),
 names_g=c("ACITR first", "APITR first"), color_g=c("red", "blue"), n.rows=3)

R.fig1$F.plot + labs(x="Quarters", color="Ordering", fill="Ordering")
R.fig2$F.plot + labs(x="Quarters", color="Ordering", fill="Ordering")

---

Impulse Response Functions for panel SVAR models

Description

Calculates impulse response functions for panel VAR objects.

Usage

## S3 method for class 'pvarx'
irf(x, ..., n.ahead = 20, normf = NULL, w = NULL, MG_IRF = TRUE)

Arguments

x	Panel VAR object of class 'pid' or 'pvarx' or a list of VAR objects that will be coerced to 'varx'.
...	Currently not used.
n.ahead	Integer. Number of periods to consider after the initial impulse, i.e. the horizon of the IRF.
normf	Function. A given function that normalizes the $K \times S$ input-matrix into an output matrix of same dimension. See the example in id.iv for the normalization of Jentsch and Lunsford (2021) that fixes the size of the impact response.
w	Numeric, logical, or character vector. $N$ numeric elements weighting the individual coefficients, or names or $N$ logical elements selecting a subset from the individuals $i = 1, \ldots, N$ for the MG estimation. If NULL (the default), all $N$ individuals are included without weights.
MG_IRF	Logical. If TRUE (the default), the mean-group of individual IRF is calculated in accordance with Gambacorta et al. (2014). If FALSE, the IRF is calculated for the mean-group of individual VAR estimates.

Value

A list of class 'svarirf' holding the impulse response functions as a 'data.frame'.

References

Luetkepohl, H. (2005): New Introduction to Multiple Time Series Analysis, Springer, 2nd ed.

Gambacorta L., Hofmann B., and Peersman G. (2014): "The Effectiveness of Unconventional Monetary Policy at the Zero Lower Bound: A Cross-Country Analysis", Journal of Money, Credit and Banking, 46, pp. 615-642.

Jentsch, C., and Lunsford, K. G. (2021): "Asymptotically Valid Bootstrap Inference for Proxy SVARs", Journal of Business and Economic Statistics, 40, pp. 1876-1891.

Examples

data("PCAP")
names_k = c("g", "k", "l", "y")  # variable names
names_i = levels(PCAP$id_i)      # country names
L.data  = sapply(names_i, FUN=function(i) 
  ts(PCAP[PCAP$id_i==i, names_k], start=1960, end=2019, frequency=1), 
  simplify=FALSE)

# estimate and identify panel SVAR #
L.vars = lapply(L.data, FUN=function(x) vars::VAR(x, p=2, type="both"))
R.pid  = pid.chol(L.vars, order_k=names_k)

# calculate and plot MG-IRF #
library("ggplot2")
R.irf = irf(R.pid, n.ahead=60)
F.irf = plot(R.irf, selection=list(2:4, 1:2))
as.pplot(F.irf=F.irf, color_g="black", n.rows=3)$F.plot + guides(color="none")

---

Impulse Response Functions

Description

Calculates impulse response functions.

Usage

## S3 method for class 'varx'
irf(x, ..., n.ahead = 20, normf = NULL)

Arguments

x	VAR object of class 'varx', 'id', or any other that will be coerced to 'varx'.
...	Currently not used.
n.ahead	Integer. Number of periods to consider after the initial impulse, i.e. the horizon of the IRF.
normf	Function. A given function that normalizes the $K \times S$ input-matrix into an output matrix of same dimension. See the example in id.iv for the normalization of Jentsch and Lunsford (2021) that fixes the size of the impact response.

Value

A list of class 'svarirf' holding the impulse response functions as a 'data.frame'.

References

Luetkepohl, H. (2005): New Introduction to Multiple Time Series Analysis, Springer, 2nd ed.

Jentsch, C., and Lunsford, K. G. (2021): "Asymptotically Valid Bootstrap Inference for Proxy SVARs", Journal of Business and Economic Statistics, 40, pp. 1876-1891.

Examples

data("PCIT")
names_k = c("APITR", "ACITR", "PITB", "CITB", "GOV", "RGDP", "DEBT")
names_l = c("m_PI", "m_CI")  # proxy names
names_s = paste0("epsilon[ ", c("PI", "CI"), " ]")  # shock names
dim_p   = 4  # lag-order

# estimate and identify proxy SVAR #
R.vars = vars::VAR(PCIT[ , names_k], p=dim_p, type="const")
R.idBL = id.iv(R.vars, iv=PCIT[-(1:dim_p), names_l], S2="MR", cov_u="OMEGA")
colnames(R.idBL$B) = names_s  # labeling

# calculate and plot normalized IRF #
R.norm = function(B) B / matrix(-diag(B), nrow(B), ncol(B), byrow=TRUE)
plot(irf(R.idBL, normf=R.norm))

---

Data set for the Monetary Demand Model

Description

The data set MDEM consists of annual observations for the nominal short-term interest rate $R$ and the logarithm of the real money aggregate $m1$ and real GDP $gdp$. It covers the period 1957 to 1996 $(T=40)$ for $N=19$ countries.

Usage

data("MDEM")

Format

A long-format data panel of class 'data.frame', where the columns id_i and id_t indicate the country and year respectively.

Source

The prepared data is sourced from OECD and IMF's International Financial Statistics of the year 1998, see the open terms of use. Employed by Carrion-i-Silvestre and Surdeanu (2011:24, Ch.6.1), it has been originally compiled and described in the unpublished appendix of Mark and Sul (2003). See the related working paper of Mark and Sul (1999, Appendix B).

References

Carrion-i-Silvestre, J. L., and Surdeanu L. (2011): "Panel Cointegration Rank Testing with Cross-Section Dependence", Studies in Nonlinear Dynamics & Econometrics, 15 (4), pp. 1-43.

Mark, N. C., and Sul, D. (1999): "A Computationally Simple Cointegration Vector Estimator for Panel Data", Working Paper, Department of Economics, Ohio State University.

Mark, N. C., and Sul, D. (2003): "Cointegration Vector Estimation by Panel DOLS and Long-Run Money Demand," Oxford Bulletin of Economics and Statistics, 65, pp. 655-680.

See Also

Other data sets: ERPT, EURO, EU_w, ICAP, MERM, PCAP, PCIT

---

Data set for the Monetary Exchange Rate Model

Description

The data set MERM consists of monthly observations for the log-ratios of prices $p$, money supply $m$, and industrial production $y$ as well as the natural logarithm of nominal exchange rates against the dollar $s$. It covers the period January 1995 to December 2007 $(T=156)$ for $N=19$ countries.

Usage

data("MERM")

Format

A long-format data panel of class 'data.frame', where the columns id_i and id_t indicate the country and month respectively.

Source

The prepared data set is directly obtainable from the journal website: doi:10.1016/j.ecosta.2016.10.001. Supplementary Raw Research Data. This is open data under the CC BY 4.0 license.

References

Oersal, D. D. K., and Arsova, A. (2017): "Meta-Analytic Panel Cointegrating Rank Tests for Dependent Panels", Econometrics and Statistics, 2, pp. 61-72.

See Also

Other data sets: ERPT, EURO, EU_w, ICAP, MDEM, PCAP, PCIT

---

Data set on Public Capital Stocks

Description

The data set PCAP consists of annual observations for

• the governmental capital stocks $G$ and their logarithm $g$,

• the private capital stocks $K$ and their logarithm $k$,

• the total hours worked $L$ and their logarithm $l$, and

• the real GDP $Y$ and its logarithm $y$.

It is constructed from the annual observations for

• the governmental investments $IG$,

• the private non-residential investments and capital stocks $IB$ and $B$,

• the private housing investments and capital stocks $IH$ and $H$, and

• the persons employed $ET$ and hours worked per person $HRS$.

It covers the period 1960 to 2019 $(T=60)$ for $N=23$ OECD countries. All monetary values are given in US-Dollars at 2005 prices, i.e. constant PPP.

Usage

data("PCAP")

Format

A long-format data panel of class 'data.frame', where the columns id_i and id_t indicate the country and year respectively.

Source

Own compilation based on data from PWT, Eurostat, and OECD's Economic Outlook. Capital stocks are derived by the Perpetual Inventory Method as described by Kamps (2006). This is open data under the CC BY 4.0 license.

References

Empting, L. F. T., and Herwartz, H. (2025): "Revisiting the 'Productivity of Public Capital': VAR Evidence on the Heterogeneous Dynamics in a New Panel of 23 OECD Countries".

Feenstra, R. C., Inklaar, R., and Timmer, M. P. (2015): "The Next Generation of the Penn World Table", American Economic Review, 105, pp. 3150-3182.

Kamps, C. (2006): "New Estimates of Government Net Capital Stocks for 22 OECD Countries, 1960-2001", IMF Staff Papers, 53, pp. 120-150.

See Also

Other data sets: ERPT, EURO, EU_w, ICAP, MDEM, MERM, PCIT

Examples

### Latex figure: public capital stocks ###
data(PCAP)
names_i = c("DEU", "FRA", "GRC", "ITA", "NLD")  # select 5 exemplary countries
idx_i  = PCAP$id_i %in% names_i
idx_pl = c(1, 3, 5, 9, 10)
pallet = c(RColorBrewer::brewer.pal(n=11, name="Spectral")[idx_pl], "#000000")
names(pallet) = c(names_i, "\\O23")
breaks = factor(c(names_i, "\\O23"), levels=c(names_i, "\\O23"))
events = data.frame(
  label   = paste0("\\quad \\textbf{ ",
            c("Oil Crisis", "Oil Crisis II", "Early 1990s", "Early 2000s", 
              "Great Recession", "European Debt Crisis", "COVID-19"), " }"), 
  t_start = c(1973, 1979, 1990, 2000, 2007, 2010, 2020),
  t_end   = c(1975, 1982, 1993, 2003, 2010, 2013, 2022))

# plot events #
library("ggplot2")
F.base <- ggplot() + 
  geom_hline(yintercept = 0, color="grey") +
  geom_rect(data=events, aes(xmin=t_start, xmax=t_end, ymin=-Inf, ymax=+Inf), 
            fill="black", alpha=0.2) +
  geom_text(data=events, aes(x=(t_start+t_end)/2, y=-Inf, label=label), 
            hjust=0, angle=90, colour='white') +
  scale_x_continuous(breaks = seq(1960, 2020, 10), limits = c(1960, 2022)) +
  theme_bw(base_size=10)

# add levels #
F.G <- F.base +
  geom_line(data=PCAP[idx_i, ], aes(x=id_t, y=G/1e+12, colour=id_i, group=id_i), 
            linewidth=2) +
  stat_summary(data=PCAP[idx_i, ], aes(x=id_t, y=G/1e+12, color="\\O23"), 
               fun=mean, geom="line", linewidth=0) +
  geom_text(data=events, aes(x=(t_start+t_end)/2, y=-Inf, label=label), 
            hjust=0, angle=90, colour='white', alpha=0.6) +
  scale_colour_manual(values=pallet, breaks=breaks) +
  labs(x=NULL, y="Trillion US-\\$ at 2005 prices", colour="Country", title=NULL) +
  guides(colour=guide_legend(override.aes = list(linewidth=2))) +
  theme(legend.position="inside", legend.position.inside=c(0.01, 0.99), 
        legend.justification = c(0, 1), 
        legend.title=element_text(size=8), legend.text=element_text(size=6), 
        legend.key.width = unit(0.35, "cm"), legend.key.height = unit(0.35, "cm"))

# add growth rates #
PCAP$gG = ave(PCAP$G, PCAP$id_i, FUN=function(x) 
  c(diff(x), NA)/x*100)  # beginning-of-the-year observations!
F.gG <- F.base +
  geom_line(data=PCAP[idx_i, ], aes(x=id_t, y=gG, colour=id_i, group=id_i), 
            linewidth=2) +
  stat_summary(data=PCAP, aes(x=id_t, y=gG, color="\\O23"), 
               fun=mean, geom="line", linewidth=2) +
  geom_text(data=events, aes(x=(t_start+t_end)/2, y=-Inf, label=label), 
            hjust=0, angle=90, colour='white', alpha=0.6) +
  scale_colour_manual(values=pallet, breaks=breaks, guide="none") +
  labs(x="Year", y="Growth in \\%", colour="Country", title=NULL)

# export to Latex #

library(tikzDevice)
textwidth = 15.5/2.54  # LaTeX textwidth from "cm" into "inch"
file_fig  = file.path(tempdir(), "Fig_G.tex")

tikz(file=file_fig, width=1.2*textwidth, height=0.8*textwidth)
 # gridExtra::grid.arrange(grobs=list(F.G, F.gG), layout_matrix=cbind(1:2))
 ggpubr::ggarrange(F.G, F.gG, ncol=1, nrow=2, align="v")
dev.off()

---

Data set on Personal and Corporate Income Tax

Description

The data set PCIT consists of quarterly observations for

• the average personal income tax rates $APITR$,

• the average corporate income tax rates $ACITR$,

• the logarithm of the personal income tax base $PITB$,

• the logarithm of the corporate income tax base $CITB$,

• the logarithm of government spending $GOV$,

• the logarithm of GDP divided by population $RGDP$, and

• the logarithm of government debt held by the public divided by the GDP deflator and population $DEBT$.

Moreover, the proxies for shocks to personal $m\_PI$ and corporate $m\_CI$ income taxes are prepended, where non-zero observations from the related narratively identified shock series $T\_PI$ resp. $T\_CI$ have been demeaned. The data set covers the period Q1 1950 to Q4 2006 $(T=228)$ for the US.

Usage

data("PCIT")

Format

A time series data set of class 'data.frame', where the column id_t indicates the quarter of the year.

Source

The prepared data set is directly obtainable from openICPSR: doi:10.3886/E116190V1. Supplementary Research Data. This is open data under the CC BY 4.0 license.

References

Mertens, K., and Ravn, M. O. (2013): "The Dynamic Effects of Personal and Corporate Income Tax Changes in the United States", American Economic Review, 103, pp. 1212-1247.

Jentsch, C., and Lunsford, K. G. (2019): "The Dynamic Effects of Personal and Corporate Income Tax Changes in the United States: Comment", American Economic Review, 109, pp. 2655-2678.

Mertens, K., and Ravn, M. O. (2019): "The Dynamic Effects of Personal and Corporate Income Tax Changes in the United States: Reply", American Economic Review, 109, pp. 2679-2691.

See Also

Other data sets: ERPT, EURO, EU_w, ICAP, MDEM, MERM, PCAP

---

Panel cointegration rank tests

Description

Performs test procedures for the rank of cointegration in a panel of VAR models. First, the chosen individual procedure is applied over all $N$ individual entities for $r_{H0}=0,\ldots,K-1$. Then, the $K \times N$ individual statistics and $p$-values are combined to panel test results on each $r_{H0}$ using all combination approaches available for the chosen procedure.

Usage

pcoint.JO(
  L.data,
  lags,
  type = c("Case1", "Case2", "Case3", "Case4"),
  t_D1 = NULL,
  t_D2 = NULL,
  n.factors = FALSE
)

pcoint.BR(
  L.data,
  lags,
  type = c("Case1", "Case2", "Case3", "Case4"),
  t_D1 = NULL,
  t_D2 = NULL,
  n.iterations = FALSE
)

pcoint.SL(L.data, lags, type = "SL_trend", t_D = NULL, n.factors = FALSE)

pcoint.CAIN(L.data, lags, type = "SL_trend", t_D = NULL)

Arguments

L.data	List of 'data.frame' objects for each individual. The variables must have the same succession $k = 1,\ldots,K$ in each individual 'data.frame'.
lags	Integer or vector of integers. Lag-order of the VAR models in levels, which is either a common $p$ for all individuals or individual-specific $p_i$ for each individual. In the vector, $p_i$ must have the same succession $i = 1,\ldots,N$ as argument L.data.
type	Character. The conventional case of the deterministic term.
t_D1	List of vectors. The activating break periods $\tau$ for the period-specific deterministic regressors in $d_{1,it}$, which are restricted to the cointegration relations. The accompanying lagged regressors are automatically included in $d_{2,it}$. The $p$-values are calculated for up to two breaks resp. three sub-samples.
t_D2	List of vectors. The activating break periods $\tau$ for the period-specific deterministic regressors in $d_{2,it}$, which are unrestricted.
n.factors	Integer. Number of common factors to be subtracted for the PANIC by Arsova and Oersal (2017, 2018). Deactivated if FALSE (the default).
n.iterations	Integer. The (maximum) number of iterations for the pooled estimation of the cointegrating vectors.
t_D	List of vectors. The activating break periods $\tau$ for the period-specific deterministic regressors in $d_{it}$ of the SL-procedure. The accompanying lagged regressors are automatically included in $d_{it}$. The $p$-values are calculated for up to two breaks resp. three sub-samples.

Value

A list of class 'pcoint' with elements:

panel	List for the panel test results, which contains one matrix for the test statistics and one for the $p$-values.
individual	List for the individual test results, which contains one matrix for the test statistics and one for the $p$-values.
CSD	List of measures for cross-sectional dependency. NULL if a first generation test has been used.
args_pcoint	List of characters and integers indicating the panel cointegration test and specifications that have been used.
beta_H0	List of matrices, which comprise the pooled cointegrating vectors for each rank $r_{H0}$. NULL if any other test than BR has been used.

Functions

• pcoint.JO(): based on the Johansen procedure.

• pcoint.BR(): based on the pooled two-step estimation of the cointegrating vectors.

• pcoint.SL(): based on the Saikkonen-Luetkepohl procedure.

• pcoint.CAIN(): accounting for correlated probits between the individual SL-procedures.

References

Larsson, R., Lyhagen, J., and Lothgren, M. (2001): "Likelihood-based Cointegration Tests in Heterogeneous Panels", Econometrics Journal, 4, pp. 109-142.

Choi, I. (2001): "Unit Root Tests for Panel Data", Journal of International Money and Finance, 20, pp. 249-272.

Arsova, A., and Oersal, D. D. K. (2018): "Likelihood-based Panel Cointegration Test in the Presence of a Linear Time Trend and Cross-Sectional Dependence", Econometric Reviews, 37, pp. 1033-1050.

Breitung, J. (2005): "A Parametric Approach to the Estimation of Cointegration Vectors in Panel Data", Econometric Reviews, 24, pp. 151-173.

Oersal, D. D. K., and Droge, B. (2014): "Panel Cointegration Testing in the Presence of a Time Trend", Computational Statistics & Data Analysis, 76, pp. 377-390.

Oersal, D. D. K., and Arsova, A. (2017): "Meta-Analytic Cointegrating Rank Tests for Dependent Panels", Econometrics and Statistics, 2, pp. 61-72.

Arsova, A., and Oersal, D. D. K. (2018): "Likelihood-based Panel Cointegration Test in the Presence of a Linear Time Trend and Cross-Sectional Dependence", Econometric Reviews, 37, pp. 1033-1050.

Hartung, J. (1999): "A Note on Combining Dependent Tests of Significance", Biometrical Journal, 41, pp. 849-855.

Arsova, A., and Oersal, D. D. K. (2021): "A Panel Cointegrating Rank Test with Structural Breaks and Cross-Sectional Dependence", Econometrics and Statistics, 17, pp. 107-129.

Examples

### reproduce Oersal,Arsova 2017:67, Ch.5 ###
data("MERM")
names_k = colnames(MERM)[-(1:2)] # variable names
names_i = levels(MERM$id_i)      # country names
L.data  = sapply(names_i, FUN=function(i) 
   ts(MERM[MERM$id_i==i, names_k], start=c(1995, 1), frequency=12), 
   simplify=FALSE)

# Oersal,Arsova 2017:67, Tab.5 #
R.lags = c(2, 2, 2, 2, 1, 2, 2, 4, 2, 3, 2, 2, 2, 2, 2, 1, 1, 2, 2)
names(R.lags) = names_i  # individual lags by AIC (lag_max=4)
n.factors = 8  # number of common factors by Onatski's (2010) criterion
R.pcsl = pcoint.SL(L.data, lags=R.lags, n.factors=n.factors, type="SL_trend")
R.pcjo = pcoint.JO(L.data, lags=R.lags, n.factors=n.factors, type="Case4")

# Oersal,Arsova 2017:67, Tab.6 #
R.Ftsl = coint.SL(y=R.pcsl$CSD$Ft, dim_p=2, type_SL="SL_trend")  # lag-order by AIC
R.Ftjo = coint.JO(y=R.pcsl$CSD$Ft, dim_p=2, type="Case4")

### reproduce Oersal,Arsova 2016:13, Ch.6 ###
data("ERPT")
names_k = c("lpm5", "lfp5", "llcusd")  # variable names for "Chemicals and related products"
names_i = levels(ERPT$id_i)[c(1,6,2,5,4,3,7)]  # ordered country names
L.data  = sapply(names_i, FUN=function(i) 
   ts(ERPT[ERPT$id_i==i, names_k], start=c(1995, 1), frequency=12), 
   simplify=FALSE)

# Oersal,Arsova 2016:21, Tab.6 (only for individual results) #
R.lags = c(3, 3, 3, 4, 3, 3, 3); names(R.lags)=names_i  # lags of VAR model by MAIC
R.cain = pcoint.CAIN(L.data, lags=R.lags, type="SL_trend")
R.pcsl = pcoint.SL(L.data,   lags=R.lags, type="SL_trend")

# Oersal,Arsova 2016:22, Tab.7/8 #
R.lags = c(3, 3, 3, 4, 4, 3, 4); names(R.lags)=names_i  # lags of VAR model by MAIC
R.t_D  = list(t_break=89)  # a level shift and trend break in 2002_May for all countries
R.cain = pcoint.CAIN(L.data, lags=R.lags, t_D=R.t_D, type="SL_trend")

---

Recursive identification of panel SVAR models via Cholesky decomposition

Description

Given an estimated panel of VAR models, this function uses the Cholesky decomposition to identify 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 decomposition of the least squares covariance matrix $\Sigma_{u,i} = B_i B_i'$.

Usage

pid.chol(x, order_k = NULL)

Arguments

x	An object of class 'pvarx' or a list of VAR objects that will be coerced to 'varx'. Estimated panel of VAR objects.
order_k	Vector. Vector of characters or integers specifying the assumed structure of the recursive causality. Change the causal ordering in the instantaneous effects without permuting variables and re-estimating the VAR model.

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.
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.

References

Luetkepohl, H. (2005): New Introduction to Multiple Time Series Analysis, Springer, 2nd ed.

Sims, C. A. (2008): "Macroeconomics and Reality", Econometrica, 48, pp. 1-48.

See Also

Other panel identification functions: pid.cvm(), pid.dc(), pid.grt(), pid.iv()

Examples

data("PCAP")
names_k = c("g", "k", "l", "y")  # variable names
names_i = levels(PCAP$id_i)      # country names
L.data  = sapply(names_i, FUN=function(i) 
  ts(PCAP[PCAP$id_i==i, names_k], start=1960, end=2019, frequency=1), 
  simplify=FALSE)

# estimate and identify panel SVAR #
L.vars = lapply(L.data, FUN=function(x) vars::VAR(x, p=2, type="both"))
R.pid  = pid.chol(L.vars, order_k=names_k)

---

Independence-based identification of panel SVAR models via Cramer-von Mises (CVM) distance

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 Cramer-von Mises distance (Genest and Remillard, 2004), determines least dependent structural shocks. The minimum is obtained by a two step optimization algorithm similar to the technique described in Herwartz and Ploedt (2016).

Usage

pid.cvm(
  x,
  combine = c("group", "pool", "indiv"),
  n.factors = NULL,
  dd = NULL,
  itermax = 500,
  steptol = 100,
  iter2 = 75
)

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.
dd	Object of class 'indepTestDist' generated by 'indepTest' from package 'copula'. Simulated independent sample(s) of the same size as the data. If the sample sizes $T_i - p_i$ differ across '$i$', the strictly separated identification requires a list of $N$ individual 'indepTestDist'-objects with respective sample sizes. If NULL (the default), a suitable object will be calculated during the call of pid.cvm.
itermax	Integer. Maximum number of iterations for DEoptim.
steptol	Numeric. Tolerance for steps without improvement for DEoptim.
iter2	Integer. Number of iterations for the second optimization.

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.

References

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

Genest, C., and Remillard, B. (2004): "Tests of Independence and Randomness based on the Empirical Copula Process", Test, 13, pp. 335-370.

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.

Herwartz, H. (2018): "Hodges Lehmann detection of structural shocks - An Analysis of macroeconomic Dynamics in the Euro Area", Oxford Bulletin of Economics and Statistics, 80, pp. 736-754.

Herwartz, H., and Ploedt, M. (2016): "The Macroeconomic Effects of Oil Price Shocks: Evidence from a Statistical Identification Approach", Journal of International Money and Finance, 61, pp. 30-44.

See Also

... the individual id.cvm by Lange et al. (2021) in svars. Note that pid.cvm relies on a modification of their procedure and thus performs ICA on the pre-whitened shocks 'eps0' directly.

Other panel identification functions: pid.chol(), pid.dc(), pid.grt(), pid.iv()

Examples

# select minimal or full example #
is_min = TRUE
idx_i  = ifelse(is_min, 1, 1:14)

# load and prepare data #
data("EURO")
data("EU_w")
names_i = names(EURO[idx_i+1])  # country names (#1 is EA-wide aggregated data)
idx_k   = 1:4   # endogenous variables in individual data matrices
idx_t   = 1:76  # periods from 2001Q1 to 2019Q4
trend2  = idx_t^2

# individual VARX models with common lag-order p=2 #
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]))
L.vars = sapply(names_i, FUN=function(i) 
  vars::VAR(L.data[[i]], p=2, type="both", exogen=L.exog[[i]]), 
  simplify=FALSE)

# identify under common orthogonal matrix (with pooled sample size (T-p)*N) #
S.pind = copula::indepTestSim(n=(76-2)*length(names_i), p=length(idx_k), N=100)
R.pcvm = pid.cvm(L.vars, dd=S.pind, combine="pool")
R.irf  = irf(R.pcvm, n.ahead=50, w=EU_w)
plot(R.irf, selection=list(1:2, 3:4))

# identify individually (with same sample size T-p for all 'i') #
S.pind = copula::indepTestSim(n=(76-2), p=length(idx_k), N=100)
R.pcvm = pid.cvm(L.vars, dd=S.pind, combine="indiv")
R.irf  = irf(R.pcvm, n.ahead=50, w=EU_w)
plot(R.irf, selection=list(1:2, 3:4))

---

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
)

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.

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.

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.

See Also

Other panel identification functions: pid.chol(), pid.cvm(), pid.grt(), pid.iv()

Examples

### replicate Herwartz,Wang 2024:630, Ch.4 ###

# 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)

---