I. Introduction

The main goal of this assignment is to replicate Figure 2 of Gali (1999) AER. To do this we will need to construct a structural VAR model with the long run restriction suggested by Blanchard and Quah. The remainder of this paper is as follows. Section II describes the data. Section III performs our methodology. Part A performs an ERS test to look for unit roots - we end up using the first differnces of the data from section II. Part B estimates a reduced form VAR model. Part C imposes the long run restrictions of Blanchard and Quah to obtain a structural VAR. Part D interprets the meaning of the SVAR obtained in part C. Part E plots the cumulative IRFs from the SVAR model in part C. And finally in part F we compare our cumulative IRFs from part E to those of Gali (1999) AER. For further examination of our SVAR we also include the FEVD in part G of Section III. Section IV concludes.

II. The Data

We download the data from FRED using the Quandl function. We then plot the time series of these two data sets to get an idea of what these data sets visually look like.

Quandl.api_key("75HFthvbcxsrteNWupus")

data.real.output <- Quandl("FRED/OPHNFB", type = 'zoo')
data.hours <- Quandl("FRED/HOANBS", type = 'zoo')

y1 <- data.real.output 
y2 <- data.hours

III. Methodology

Part A - ERS - Checking for a Unit Root

As shown in Section II, the time series of our data might be non-linear, therefore we begin by transforming the data using two different transformation techniques: (1) Log transformation and (2) First differencing. We then conduct an ERS test for both transformations and for both time series to check for the precense of a unit root.

The null hypothesis of this test is that the time series has a unit root. A unit root implies that the time series will be difference stationary, that is, the mean is non-constant, non-zero, and is stochastic. The alternative hypothesis is that the time series is stationary, level stationary, or trend stationary depending upon which specification you use in the test. A stationary time series under this test represents a time series with a constant mean of zero. A level stationary time series is a time series with a non-zero but constant mean. A trend stationary time series has a mean that is non-constant, and non-zero, but is deterministic. In short, this test assumes the time series is difference stationary and tests to see if there is enough evidence to reject this hypothesis and conclude that the time series is trend stationary. The code we use to run these tests is given below and the results follow.

log.y1 <- 100*log(y1)
log.y2 <- 100*log(y2)

first.diff.y1 <- diff(y1, differences = 1)
first.diff.y2 <- diff(y2, differences = 1)

ERS.1 <- summary(ur.ers(log.y1, type = "DF-GLS", model = "trend"))
ERS.2 <- summary(ur.ers(log.y2, type = "DF-GLS", model = "trend"))

ERS.3 <- summary(ur.ers(first.diff.y1, type = "DF-GLS", model = "trend"))
ERS.4 <- summary(ur.ers(first.diff.y2, type = "DF-GLS", model = "trend"))
ERS.1       # Fail to reject null

############################################### 
# Elliot, Rothenberg and Stock Unit Root Test # 
############################################### 

Test of type DF-GLS 
detrending of series with intercept and trend 


Call:
lm(formula = dfgls.form, data = data.dfgls)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.99346 -0.51880 -0.09269  0.45875  2.75287 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
yd.lag       -0.006828   0.008578  -0.796    0.427
yd.diff.lag1  0.032803   0.061000   0.538    0.591
yd.diff.lag2  0.085846   0.059267   1.448    0.149
yd.diff.lag3 -0.021792   0.057598  -0.378    0.705
yd.diff.lag4 -0.025814   0.057341  -0.450    0.653

Residual standard error: 0.8382 on 270 degrees of freedom
Multiple R-squared:  0.01157,   Adjusted R-squared:  -0.006738 
F-statistic: 0.6319 on 5 and 270 DF,  p-value: 0.6756


Value of test-statistic is: -0.796 

Critical values of DF-GLS are:
                 1pct  5pct 10pct
critical values -3.48 -2.89 -2.57
ERS.2       # Fail to reject null

############################################### 
# Elliot, Rothenberg and Stock Unit Root Test # 
############################################### 

Test of type DF-GLS 
detrending of series with intercept and trend 


Call:
lm(formula = dfgls.form, data = data.dfgls)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.31024 -0.35982 -0.01214  0.41156  2.63951 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
yd.lag       -0.014995   0.007663  -1.957   0.0514 .  
yd.diff.lag1  0.662323   0.059896  11.058   <2e-16 ***
yd.diff.lag2 -0.000947   0.072072  -0.013   0.9895    
yd.diff.lag3 -0.009683   0.072056  -0.134   0.8932    
yd.diff.lag4 -0.113993   0.060497  -1.884   0.0606 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6814 on 270 degrees of freedom
Multiple R-squared:  0.4322,    Adjusted R-squared:  0.4217 
F-statistic:  41.1 on 5 and 270 DF,  p-value: < 2.2e-16


Value of test-statistic is: -1.9567 

Critical values of DF-GLS are:
                 1pct  5pct 10pct
critical values -3.48 -2.89 -2.57
ERS.3       # Rejects the null at the 1% significance level

############################################### 
# Elliot, Rothenberg and Stock Unit Root Test # 
############################################### 

Test of type DF-GLS 
detrending of series with intercept and trend 


Call:
lm(formula = dfgls.form, data = data.dfgls)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.38068 -0.34340 -0.07534  0.21912  1.63935 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
yd.lag       -0.62856    0.10547  -5.960  7.9e-09 ***
yd.diff.lag1 -0.27088    0.10084  -2.686  0.00767 ** 
yd.diff.lag2 -0.12760    0.09387  -1.359  0.17522    
yd.diff.lag3 -0.06366    0.08196  -0.777  0.43798    
yd.diff.lag4 -0.01179    0.06054  -0.195  0.84568    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4971 on 269 degrees of freedom
Multiple R-squared:  0.4506,    Adjusted R-squared:  0.4404 
F-statistic: 44.13 on 5 and 269 DF,  p-value: < 2.2e-16


Value of test-statistic is: -5.9596 

Critical values of DF-GLS are:
                 1pct  5pct 10pct
critical values -3.48 -2.89 -2.57
ERS.4       # Rejects the null at the 1% significance level 

############################################### 
# Elliot, Rothenberg and Stock Unit Root Test # 
############################################### 

Test of type DF-GLS 
detrending of series with intercept and trend 


Call:
lm(formula = dfgls.form, data = data.dfgls)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.01209 -0.30345  0.02839  0.34190  1.91103 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
yd.lag       -0.363646   0.062214  -5.845 1.46e-08 ***
yd.diff.lag1 -0.001562   0.068460  -0.023   0.9818    
yd.diff.lag2  0.124534   0.065844   1.891   0.0597 .  
yd.diff.lag3  0.055696   0.064604   0.862   0.3894    
yd.diff.lag4 -0.065018   0.060827  -1.069   0.2861    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5023 on 269 degrees of freedom
Multiple R-squared:  0.1897,    Adjusted R-squared:  0.1746 
F-statistic: 12.59 on 5 and 269 DF,  p-value: 5.331e-11


Value of test-statistic is: -5.8451 

Critical values of DF-GLS are:
                 1pct  5pct 10pct
critical values -3.48 -2.89 -2.57

The first two tests, examine the ERS test statistics of the log transformed data. Both tests fail to reject the null hypothesis which implies that the log-transformed data contains a unit root and therefore is not approximately weakly stationary. The last two tests examine the ERS test statistic of the first difference tranformed data. Both tests reject the null hypothesis at the 1% level, which gives us ground to assume that the first difference of both of our data sets are approximately weakly stationary.

Due to the results of the ERS tests, for the remainder of this assignment we will be using the first differnces of our data sets. Hence, we redefine our data as below for convience of names.

real.output <- first.diff.y1
hours <- first.diff.y2

Part B - Estimating the Bivariate Reduced Form VAR Model

We now estimate a bivariate reduced form vector autoregression model (VAR) for \(\mathbf{y_t} = (\Delta y_{1,t}, \Delta y_{2,t})'\) for the sample period (1947Q2 - 2016Q4). (Note: we lose one quarter of data, 1947Q1 because we took the first difference). We create the vector \(\mathbf{y_t}\) simply by combining \(y_{1,t}\) and \(y_{2,t}\) as columns over each quarter in time across the sample. We also use the ic input to specify that we want the number of lags in the VAR model to be chosen based upon the AIC information criteria statistic. We output the results in Table I which show that the AIC statistic chose two lags for our VAR system.

y <- cbind(real.output, hours)
var.p <- VAR(y, ic = "AIC", lag.max = 12, type = "const")
Table I: Reduced form Bivariate VAR Model
Dependent variable:
real.output hours
real.output.l1 -0.014 0.181***
(0.058) (0.062)
hours.l1 0.011 0.611***
(0.056) (0.060)
real.output.l2 0.071 0.083
(0.058) (0.062)
hours.l2 -0.232*** 0.081
(0.057) (0.061)
const 0.334*** -0.006
(0.039) (0.042)
Observations 277 277
R2 0.103 0.452
Adjusted R2 0.089 0.444
Residual Std. Error (df = 272) 0.463 0.496
F Statistic (df = 4; 272) 7.776*** 56.038***
Note: p<0.1; p<0.05; p<0.01

Part C - Structural VAR Model

In this section we want to study the effects of technology shocks and demand shocks on hours worked. To do so we use a Blanchard and Quah approach to obtain an structural vector autoregressive model (SVAR) model where we impose the condition that demand shocks do not affect real output per hour, \(y_{1,t}\), in the long run.

In general, the SVAR model uses economic theory to sort out contemporaneous links among the variables. SVAR models require identifying assumptions that establish causal links among varaibles, that is, models the contemoraneous interdependence between the left hand side variables (the non-lagged variables, the dependent variables in our regressions of part B). When using the reduced form VAR we are using Choleski’s decomposition, however, with an SVAR model our identifying assumptions use economic theory to impose restrictions to achieve identification of parameters of the structural VAR.

There are three types of restrictions: (1) Short run restrictions, (2) Long run restrictions, and (3) Sign restrictions. For this assignment we want to impose a long run restriction, namely that demand shocks do not affect real output per hour. We will use Blanchard and Quah’s long run restriction approach. Before we can apply Blanchard and Quah’s approach we need to de-mean the data. Below we show a plot of the time series before we de-mean the data.

Since the data does not have a mean of zero this gives us ground for us to perform a de-meaning function on the data, using the code below. A plot of the de-meaned data follows.

y <- sweep(y, 2, apply(y, 2, mean))

Now that the data is de-meaned we can re-run our reduced form VAR. The results of this model are reported in Table II.

myVAR <- VAR(y, ic = "AIC", lag.max = 12)
Table II: Reduced form and De-Meaned Bivariate VAR Model
Dependent variable:
real.output hours
real.output.l1 -0.014 0.181***
(0.058) (0.062)
hours.l1 0.011 0.611***
(0.056) (0.060)
real.output.l2 0.071 0.083
(0.058) (0.062)
hours.l2 -0.232*** 0.081
(0.057) (0.061)
const 0.003 0.001
(0.028) (0.030)
Observations 277 277
R2 0.103 0.452
Adjusted R2 0.089 0.444
Residual Std. Error (df = 272) 0.463 0.496
F Statistic (df = 4; 272) 7.776*** 56.038***
Note: p<0.1; p<0.05; p<0.01

Finally, using the BQ function we can transform our reduced form VAR into an SVAR system.

mySVAR <- BQ(myVAR)

Note: It is important to note that the actual coefficeints of the model have not changed, but that we have simply imposed a long run restriction on the demand shocks to real output per hour. These changes will be discuessed in more detail in part D.

Part D - SVAR Interpretation

We now report and interpret the contemporaneous and long run impact matrices of the SVAR model reported below.

mySVAR

SVAR Estimation Results:
======================== 


Estimated contemporaneous impact matrix:
            real.output  hours
real.output      0.3623 0.2881
hours           -0.2893 0.4028

Estimated identified long run impact matrix:
            real.output hours
real.output      0.5029 0.000
hours           -0.5074 1.306

The contemporaneous impact matrix (the first matrix in the output) reports the immediate effect of a technology and demand shock on real output per hour and the amount of hours worked upon impact. In other words, the contemporanoues impact matrix tells us the affect shocks have on our VAR variables at the exact moment when the shock occurs. The long run impact matrix (the second matrix in the output) reports the long run effect of a technology and demand shock on real output per hour and the amount of hours worked across time.

Specifically, on impact, a positive one standard deviation technology shock increases real output per hour by 0.36%. This makes sense because if technology increaes we would expect workers to be more productive (since real output per hour is a proxy for labor productivity). Similarly, a positive one standard deviation technology shock decreases the amount of hours worked by 0.29% upon impact. This also is intuitive because more technology should reduce the amount of labor needed since the technology is now completing the job that labor used to be doing. Hence, the amount of labor hours should decrease. A positive one standard deviation demand shock increases real output per hour by 0.29% and increases hours by 0.40% on impact.

Now looking at the long run matrix we see that a positive one standard deviation technology shock increases real output by 0.50% in the long run as compared to 0.36% on impact. Hence, over time we say that a positive technology shock will increase labor productivity in the long run by 0.14% more than what the shock incurred at impact. This makes sense because economies of scale and labor specialization can occur in the long run with regards to the new technology, which cannot be achieved on impact. Similarly a positive one standard deviation technology shock decreases labor hours in the long run by 0.51% which is 0.23% larger in magnitude than the on impact shock. This also makes sense because as laborers become more efficeint as using the technology less and less people will be needed to help perform the task that the technology is now performing.

However, in the long run impact matrix, the demand shock does not affect the real output per hour in the long run, as indicated by the impact matrix having a value of zero. This restriction is the imposed restriction by the Blanchard and Quah approach. A positive one standard deviation demand shock positively affects the amount of hours worked in the labor market by 1.3%, more than triple the on impact affect.

Part E - Cumulative Impulse Response Functions

In this section we plot and interpret the cumulative impulse response functions based on the structural VAR model that we made in part C. Recall that the goal of the impulse response function is to track the response of one of our variables to a one time shock. In other words, we want to know how our dynamic system (the system of our variables as endogenoues functions of one another) responds to some external change which we call an impulse (sometimes called a brief input signal). The idea is we want to know how \(y_{1,t}\) and \(y_{2,t}\) respond to some exogenous shock. Tthe impulse response analysis creates an impulse response function that measures the reaction of our endogenoues macroeconomic variables to the exogenous shock at the time of the shock and over subsequent points in time. Mathematically we say that we are obtaining the response of \(y_i\) across time to a one time increase in \(\varepsilon_{j,t}\).

Similarly, the cumulative impulse response function plots the accumulation of the impact of the shock to our variables across time, instead of the shock’s impact at a single point in time. Below is our cumulative response function, we then plot the IRF using a function called plotIRF which we include in the appendix. We compute the IRF up to 12 periods ahead, to follow Gali (1999), which we will compare in part F.

myIRF.c <- irf(mySVAR, n.ahead = 12, ci = .9, cumulative = TRUE)
par(mfrow=c(2,2), cex=.6, mar=c(4,4,2,1))
source("C:/Users/rmerrima/Dropbox/Spring 2017/ECO 5316 - Time Series Econometrics/Homework 6/plotIRF.R")
plotIRF(myIRF.c, lwd=2.5, ask=FALSE, 
vlabels=c("Real Output","Hours"), 
slabels=c("technology shock","demand shock"))

The top left graph shows that real output is positively affected by a one standard deviation technology shock upon impact and that as time goes on real output continues to be positively affected by this one time shock. Similarly, the bottom left graph shows that hours are negatively affected by a one standard deviation technology shock upon impact, however the long run impact is uncertain. Hours could be unaffected by the technology shock in the long run or could be worse (lower) over time.

The top right graph explicitly shows the imposed restrictrion that the demand shock to real output per hour is zero in the long run. The initial demand shock positively affects real output per hour but this benefit marginally decreases over time until eventually in the long run, the demand shock eventually has no net gain or loss to real output per hour. The bottom right graph shows us that the amount of hours works increases do to the initial shock and that over time hours continually increase until they level off about 8 quarters after the demand shock.

Part F - Compare to Gali (1999) AER

We now compare our cumulative impulse response function to the one constructed by Gali (1999) as published in the American Economic Review. Specifically, we compare our cumulative IRF to Figure 2 of Gali’s paper Technology, Employment, and the Business Cycle: Do Technology Shocks Explain Aggregate Fluctuations? Below is Figure 2 of his paper.

One thing to point out is that Gali used three variables in his SVAR system where we only use two. Namely, the cumulative impulse response functions that we are trying to replicateis the top row of his figure 2 and the bottom row of his figure 2. Our results are relatively similar to those of Gali (1999) especially in the pattern of the graphs. However, the magniture of our impulse response functions appear to be quite different than Gali (1999). We expect that if we were to include the GDP variable that Gali includes as our third VAR variable, then our results will almost exactly replicate his.

Part G - Forecast Error Variance Decomposition

In this section wer compute the forecast error variance decomposition of our structural VAR model from part C. Below is the code to calculate the FEVD, we then present the results in Table III and plot these results in the figure below. We follow it with an interpretation of the FEVD.

myFEVD <- fevd(mySVAR, n.ahead = 12)
Table 3: FEVD of the VAR Model
y1.e -> y1 y2.e -> y1 y1.e -> y2 y2.e -> y2
Q1 0.6127 0.3873 0.3403 0.6597
Q2 0.6128 0.3872 0.2766 0.7234
Q3 0.6140 0.3860 0.2449 0.7551
Q4 0.6036 0.3964 0.2326 0.7674
Q5 0.5958 0.4042 0.2279 0.7721
Q6 0.5920 0.4080 0.2263 0.7737
Q7 0.5905 0.4095 0.2258 0.7742
Q8 0.5899 0.4101 0.2257 0.7743
Q9 0.5897 0.4103 0.2257 0.7743
Q10 0.5897 0.4103 0.2257 0.7743
Q11 0.5897 0.4103 0.2257 0.7743
Q12 0.5896 0.4104 0.2257 0.7743 
plot(myFEVD, addbars = 4, col = c("red3", "royalblue3"))

The FEVD shows us that the first difference change in real output per hour is explained by about 60% from the technology shock and about 40% from the demand shock for both the short run and the long run. Hence, the explanation of the variation in real ouput per hour is not affected by the long run or short run affects of each shock.

Furthermore, the FEVD shows us that the first difference change in hours is explained by about 34% from the technology shock and about 66% from the demand shock in the short run. However, as we approach the long run, the explanation of the variation decreases for the technology shock and increases for the demand shock. Specifically, in the long run the technology shock now only explains about 23% of the variation in hours whereas, the demand shock explains about 77% of the variation in hours.

IV. Conclusion

In this paper we have successfully replicated a portion of the Gali (1999) AER paper. We have also shown the affects of technology and demand shocks on real output per hour and the amount of labor hours worked in the United States. We have also succesfully constructed a reduced form and structural VAR model to endogenously determine these two variables across time.

Appendix - plotIRF Function

Below is the code of the modified plotIRF function used in this code. See the original at: vars plotIRF

# Sligthly modified version of plot.irf function from vars package:
# https://cran.r-project.org/web/packages/vars/

plotIRF <-
  function (x, vnames = NULL, snames = NULL, vlabels = NULL, slabels = NULL,
            main = NULL, sub = NULL, lty = NULL, lwd = NULL, col = NULL, ylim = NULL,
            ylab = NULL, xlab = NULL, nc, mar.multi = c(0, 4, 0, 4),
            oma.multi = c(6, 4, 6, 4), adj.mtext = NA, padj.mtext = NA, col.mtext = NA, ...)
  {
    # op <- par(no.readonly = TRUE)
    # on.exit(par(op))
    ##
    ## Check arguments
    ##
    inames <- x$impulse
    if (is.null(snames)) {
      snames <- inames
    }
    else {
      snames <- as.character(snames)
      if (!(all(snames %in% inames))) { stop("\nInvalid shock name(s) supplied.\n") }
      else { inames <- snames }
    }
    nvi <- length(inames)
    rnames <- x$response
    if (is.null(vnames)) {
      vnames <- rnames
    }
    else {
      vnames <- as.character(vnames)
      if (!(all(vnames %in% rnames))) { stop("\nInvalid variable name(s) supplied.\n") }
      else { rnames <- vnames }
    }
    nvr <- length(rnames)
    if (is.null(slabels)) {
      slabels <- snames
    }
    else {
      if (!(length(slabels) == length(snames))) { 
        stop("\nNUmber of labels for shocks in slabels does not math number of shocks in snames.\n") }
    }
    if (is.null(vlabels)) {
      vlabels <- vnames
    }
    else {
      if (!(length(slabels) == length(snames))) { 
        stop("\nNUmber of labels for shocks in slabels does not math number of shocks in snames.\n") }
    }
    
    ##
    ## Presetting certain plot-argument
    ##
    ifelse(is.null(lty), lty <- c(1, 1, 2, 2), lty <- rep(lty, 4)[1:4])
    ifelse(is.null(lwd), lwd <- c(1, 1, 1, 1), lwd <- rep(lwd, 4)[1:4])
    ifelse(is.null(col), col <- c("black", "gray", "red", "red"), col <- rep(col, 4)[1:4])
    
    ##
    ## Extract data for variable rname from varirf object
    ##
    dataplot.r <- function(x, rname){
      impulses <- NULL
      for(j in 1:length(x$irf)){
        impulses <- cbind(impulses, x$irf[[j]][,rname])
        colnames(impulses)[j] <- names(x$irf)[j]
      }
      range <- range(impulses)
      upper <- NULL
      lower <- NULL
      if(x$boot){
        for(j in 1:length(x$irf)){
          upper <- cbind(upper, x$Upper[[j]][,rname])
          lower <- cbind(lower, x$Lower[[j]][,rname])
          colnames(upper)[j] <- names(x$irf)[j]
          colnames(lower)[j] <- names(x$irf)[j]
        }
        range <- range(cbind(impulses, upper, lower))
      }
      if ((x$model == "varest") || (x$model == "vec2var")) {
        text1 <- "Impulse Response"
      } else if (x$model == "svarest") {
        text1 <- "SVAR Impulse Response"
      } else if (x$model == "svecest") {
        text1 <- "SVECM Impulse Response"
      }
      if (x$cumulative)  text1 <- paste("Cumulative", text1, sep = " ")
      text2 <- ""
      if (x$boot) text2 <- paste((1 - x$ci) * 100, "% Bootstrap CI, ", x$runs, "runs")
      
      result <- list(impulses = impulses, upper = upper, 
                     lower = lower, range = range, 
                     text1 = text1, text2 = text2)
      return(result)
    }
    ##
    ## Plot function for irf per impulse and response
    ##
    plot.single <- function(x, iname, rname, ylabel, slabel,...) {
      x$text1 <- paste(x$text1, "from", slabel, sep = " ")
      ifelse(is.null(main), main <- x$text1, main <- main)
      ifelse(is.null(sub), sub <- x$text2, sub <- sub)
      xy <- xy.coords(x$impulses[, iname])
      ifelse(is.null(xlab), xlabel <- "", xlabel <- xlab)
      ifelse(is.null(ylim), ylim <- x$range, ylim <- ylim)
      plot(xy, type = "l", ylim = ylim, 
           col = col[1], lty = lty[1], 
           lwd = lwd[1], axes = FALSE, 
           ylab = paste(ylabel), xlab = paste(xlab), ...)
      title(main = main, sub = sub, ...)
      axis(1, at = xy$x, labels = c(0:(length(xy$x) - 1)))
      axis(2, ...)
      box()
      if (!is.null(x$upper)) lines(x$upper[, iname], col = col[3], lty = lty[3], lwd = lwd[3])
      if (!is.null(x$lower)) lines(x$lower[, iname], col = col[3], lty = lty[3], lwd = lwd[3])
      abline(h = 0, col = col[2], lty = lty[2], lwd = lwd[2])
    }
    ##
    ## Plot IRFs
    ##
    for(i in 1:nvr){
      dp.r <- dataplot.r(x, rname = rnames[i])
      for(j in 1:nvi){
        plot.single(dp.r, iname = inames[j], rname = rnames[i], 
                    ylabel = vlabels[i], slabel = slabels[j], ...)
        # if (nvr > 1) par(ask = TRUE)
      }
    }
  }