Software 14.7.1 (Indirect least squares and 2SLS in a market). A simulated market with the structure of Example 14.1.2, true demand slope \(-1.5\) and supply slope \(2\); the function tsls is that of Software 13.5.1:

> set.seed(2); n <- 500
> inc <- rnorm(n, 10, 4); w <- rnorm(n, 5, 2)  # demand and supply shifters
> ud <- rnorm(n, 0, 3); us <- rnorm(n, 0, 3)
> p <- (40 + 0.8 * inc + 3 * w + ud - us) / 3.5 # market-clearing price
> q <- 50 - 1.5 * p + 0.8 * inc + ud            # on the demand curve
> rf.p <- coef(lm(p ~ inc + w)); rf.q <- coef(lm(q ~ inc + w))
> rbind(p = rf.p, q = rf.q)                     # the reduced form
  (Intercept)       inc          w
p    11.51861 0.2291435  0.8411853
q    32.58003 0.4768458 -1.2580431
> c(ILS.demand = rf.q[["w"]] / rf.p[["w"]],
+   ILS.supply = rf.q[["inc"]] / rf.p[["inc"]])
ILS.demand ILS.supply
 -1.495560   2.080992
> tsls <- function(y, X, Z) {
+   Xh <- Z %*% solve(crossprod(Z), crossprod(Z, X))
+   b <- solve(crossprod(Xh), crossprod(Xh, y))
+   u <- y - X %*% b; s2 <- sum(u^2) / (length(y) - ncol(X))
+   list(b = drop(b), se = sqrt(diag(s2 * solve(crossprod(Xh)))))
+ }
> Z <- cbind(1, inc, w)
> dem <- tsls(q, cbind(1, p, inc), Z); sup <- tsls(q, cbind(1, p, w), Z)
> round(rbind(demand = dem$b, supply = sup$b), 3)
                   p    inc
demand 49.807 -1.496  0.820
supply  8.610  2.081 -3.009
> round(c(OLS.demand = coef(lm(q ~ p + inc))[["p"]],
+         OLS.supply = coef(lm(q ~ p + w))[["p"]]), 3)
OLS.demand OLS.supply
    -0.926      0.882

Reading the output. The ratio of the reduced-form coefficients of \(w\) gives the demand slope, \(-1.258/0.841 = -1.496\), and the ratio of the income coefficients gives the supply slope, \(0.477/0.229 = 2.081\). 2SLS on each equation, with all the exogenous variables as instruments, reproduces these slopes exactly, as it must for exactly identified equations. Least squares gives \(-0.93\) and \(0.88\): each structural regression mixes the two curves.

Stata: reg3 (q p inc) (q p w), 2sls for 2SLS on each equation, or reg3 alone for 3SLS. EViews: create a System object with the two equations and the instrument list, then Estimate by 2SLS or 3SLS.

Software 14.7.2 (SUR: front and rear seat passengers). The Seatbelts data also record monthly front-seat and rear-seat passengers killed or seriously injured. The 1983 law required front-seat occupants only to wear belts. Monthly shocks such as weather affect both groups, so their equations, with the regressors of Software 8.7.1 plus the trend, form a natural SUR system. The block-diagonal stacked matrix is built with Matrix::bdiag:

> sb <- data.frame(Seatbelts); sb$t <- 1:192
> sb$month <- factor(cycle(Seatbelts[, "drivers"]))
> X1 <- model.matrix(~ log(kms) + log(PetrolPrice) + law + month + t, sb)
> X2 <- X1[, colnames(X1) != "law"]             # rear equation without law
> y1 <- log(sb$front); y2 <- log(sb$rear)
> sur <- function(ys, Xs) {                     # feasible GLS on the stack
+   n <- length(ys[[1]]); X <- as.matrix(Matrix::bdiag(Xs)); y <- unlist(ys)
+   E <- sapply(1:2, function(i) lm.fit(Xs[[i]], ys[[i]])$residuals)
+   W <- kronecker(solve(crossprod(E) / n), diag(n))
+   V <- solve(t(X) %*% W %*% X); b <- drop(V %*% t(X) %*% W %*% y)
+   list(b = b, se = sqrt(diag(V)), r = cor(E)[1, 2])
+ }
> j <- which(colnames(X1) == "law")             # position of law in eq. 1
> ols1 <- coef(summary(lm(y1 ~ X1 - 1)))[j, 1:2]
> ols2 <- coef(summary(lm(y2 ~ X1 - 1)))[j, 1:2]
> s1 <- sur(list(y1, y2), list(X1, X1))         # same regressors
> s2 <- sur(list(y1, y2), list(X1, X2))         # law excluded from rear
> round(rbind(OLS.front = ols1, OLS.rear = ols2,
+   SUR.same = c(s1$b[j], s1$se[j]), SUR.restr = c(s2$b[j], s2$se[j])), 4)
          Estimate Std. Error
OLS.front  -0.2228     0.0248
OLS.rear    0.1128     0.0275
SUR.same   -0.2228     0.0238
SUR.restr  -0.2886     0.0175
> c(r = s1$r, LM.diag = 192 * s1$r^2)           # Breusch-Pagan, diagonal
         r    LM.diag
 0.7079925 96.2406404

Reading the output. The residual correlation between the two equations is \(0.71\), and the Breusch–Pagan statistic \(96.2\) on \(1\) degree of freedom rejects a diagonal \(\Sigma \) overwhelmingly. With the law dummy in both equations the regressors are the same, and SUR reproduces least squares exactly, \(-0.2228\), as Theorem 14.4.1(ii) says. (The standard error is smaller only because SUR divides by \(n\) rather than \(n - k\): \(0.0248\sqrt {176/192} = 0.0238\).) The rear-seat equation shows a surprise: rear casualties rose by \(100(e^{0.113} - 1) = 12\%\) after the law, an increase that Harvey and Durbin (1986) also found in their analysis of the law for the Department of Transport; one explanation offered is risk compensation, drivers who feel safer driving less carefully. If one imposes the plausible-looking restriction that the law did not affect rear passengers, by dropping it from that equation, the regressors differ and SUR now “gains”: the front-seat effect becomes \(-0.289\) with a standard error of \(0.018\). But the restriction is false, and the gain is illusory. SUR has transferred the misspecification of the rear equation into the front equation, as system methods do. The rear equation itself, whose regressors are a subset of the front’s, is unchanged from least squares, as the last part of the theorem says.

Stata: sureg (lfront lkms lpetrol law i.month t) (lrear lkms lpetrol i.month t), with , corr for the Breusch–Pagan test. EViews: a System object estimated by Seemingly Unrelated Regression.

Software 14.7.3 (Lag selection and Johansen’s test). Return to the sales and leading indicator series of Software 12.6.3. Choose the order of the VAR in levels by BIC on a common sample, and compute Johansen’s statistics by the steps of Theorem 14.6.2:

> Y <- cbind(sales = as.numeric(BJsales), lead = as.numeric(BJsales.lead))
> bic <- function(k, P = 9) {                   # levels VAR(k), common sample
+   idx <- (P + 1):nrow(Y)
+   X <- cbind(1, do.call(cbind, lapply(1:k, function(j) Y[idx - j, ])))
+   E <- Y[idx, ] - X %*% solve(crossprod(X), crossprod(X, Y[idx, ]))
+   m <- length(idx); log(det(crossprod(E) / m)) + log(m) * 2 * ncol(X) / m
+ }
> round(sapply(1:8, bic), 3)
[1] -1.831 -1.928 -4.610 -4.507 -4.571 -4.575 -4.494 -4.460
> johansen <- function(Y, k) {                  # VAR(k) in levels, constant
+   D <- diff(Y); idx <- k:nrow(D); m <- length(idx)
+   Z <- cbind(1, do.call(cbind,
+                         lapply(seq_len(k - 1), function(j) D[idx - j, ])))
+   res <- function(M) M - Z %*% solve(crossprod(Z), crossprod(Z, M))
+   R0 <- res(D[idx, ]); R1 <- res(Y[idx, ])    # Y[idx, ] is Y(t-1)
+   S00 <- crossprod(R0) / m; S11 <- crossprod(R1) / m
+   S01 <- crossprod(R0, R1) / m
+   e <- eigen(solve(S11, t(S01)) %*% solve(S00, S01))
+   lam <- Re(e$values); b <- Re(e$vectors[, 1])
+   list(lambda = lam, trace = rev(cumsum(rev(-m * log(1 - lam)))),
+        beta = b / b[1])
+ }
> jo <- johansen(Y, 3); lapply(jo, round, 3)
$lambda
[1] 0.944 0.006

$trace
[1] 424.802   0.837

$beta
[1]   1.000 -18.169

> round(sapply(2:8, function(k) johansen(Y, k)$trace[1]), 1)
[1]  55.6 424.8  75.3  25.3  12.9   8.7   6.3

Reading the output. BIC chooses \(k = 3\). The trace statistic for \(r = 0\) is \(424.8\), far above \(15.41\), and for \(r \leq 1\) it is \(0.84\), below \(3.76\): there is exactly one cointegrating vector, \(\text {sales} - 18.17\,\text {lead}\), close to the Engle–Granger estimate \(17.34\) of Software 12.6.4 (which used the indicator three periods earlier). The first canonical correlation is enormous, \(\hat \lambda _1 = 0.944\), because the indicator predicts the change in sales so well. The last line is a warning: the trace statistic for \(r = 0\) falls from \(424.8\) to \(6.3\) as the lag order rises from \(3\) to \(8\), and with \(k \geq 6\) the test no longer rejects. Over-fitting the short-run dynamics absorbs the long-run relation into the lags, and the power of the test collapses.

Stata: varsoc sales lead for lag selection, vecrank sales lead, lags(3) for the trace test. EViews: in a Group, View \(\to \) Cointegration Test \(\to \) Johansen System Cointegration Test.

Software 14.7.4 (The VECM: adjustment and Granger causality). Estimate the VECM with the Johansen cointegrating vector and \(k - 1 = 2\) lagged differences, and test Granger causality in each direction (lags of the other variable together with the error correction term, which contains both):

> options(show.signif.stars = FALSE)
> D <- diff(Y); idx <- 3:nrow(D)                # VECM, 2 lagged differences
> ec <- drop(Y[idx, ] %*% jo$beta)              # sales - 18.17 lead, lagged
> X <- cbind(ec, D[idx - 1, ], D[idx - 2, ])
> colnames(X) <- c("ec", "ds1", "dl1", "ds2", "dl2")
> ds <- D[idx, 1]; dl <- D[idx, 2]
> vs <- lm(ds ~ X); vl <- lm(dl ~ X)
> round(rbind(sales = coef(summary(vs))["Xec", ],
+             lead = coef(summary(vl))["Xec", ]), 4)
      Estimate Std. Error  t value Pr(>|t|)
sales  -0.2597     0.0053 -48.7421   0.0000
lead    0.0047     0.0052   0.8973   0.3711
> gc <- function(y, keep) {                     # F test of the dropped terms
+   anova(lm(y ~ X[, keep]), lm(y ~ X))[2, c("F", "Pr(>F)")]
+ }
> rbind(lead.to.sales = gc(ds, c("ds1", "ds2")), # drops ec, dl1, dl2
+       sales.to.lead = gc(dl, c("dl1", "dl2"))) # drops ec, ds1, ds2
                      F Pr(>F)
lead.to.sales 1014.8068 <2e-16
sales.to.lead    1.1457 0.3329

Reading the output. Sales adjust strongly to the disequilibrium: \(\hat \alpha = -0.260\) with \(t = -48.7\), so a quarter of any gap between sales and \(18.17\) times the indicator closes each period. The indicator does not adjust at all (\(\hat \alpha = 0.005\), \(t = 0.90\)): it is weakly exogenous, and it carries the common stochastic trend that sales follow. Granger causality is overwhelming from the indicator to sales (\(F = 1015\)) and absent in the other direction (\(p = 0.33\)). The indicator deserves its name.

Stata: vec sales lead, lags(3), then vargranger. EViews: estimate a VAR object as a Vector Error Correction model; View \(\to \) Lag Structure \(\to \) Granger Causality.

Software 14.7.5 (Impulse responses). Responses of both variables, in levels, to an orthogonalised shock to the indicator, computed by the recursion of Theorem 14.5.2(iii) from the VAR(3) in levels. The indicator is ordered first, so its shock may move sales in the same period:

> k <- 3; idx <- (k + 1):nrow(Y)                # levels VAR(3)
> XV <- cbind(1, do.call(cbind, lapply(1:k, function(j) Y[idx - j, ])))
> B <- solve(crossprod(XV), crossprod(XV, Y[idx, ]))
> E <- Y[idx, ] - XV %*% B; S <- crossprod(E) / length(idx)
> A <- lapply(1:k, function(j) t(B[2 * j + 0:1, ]))      # A1, A2, A3
> Psi <- list(diag(2))                          # MA coefficient matrices
> for (h in 1:15) Psi[[h + 1]] <- Reduce(`+`,
+   lapply(1:min(h, k), function(j) A[[j]] %*% Psi[[h + 1 - j]]))
> L <- t(chol(S[2:1, 2:1]))[2:1, 2:1]           # Cholesky, lead ordered first
> resp <- sapply(Psi, function(M) M %*% L[, 2]) # responses to a lead shock
> dimnames(resp) <- list(c("sales", "lead"), paste0("h", 0:15))
> round(resp[, c(1:7, 11, 16)], 2)
        h0   h1   h2   h3   h4   h5   h6  h10  h15
sales 0.01 0.01 0.03 1.31 1.57 1.83 2.02 2.67 2.87
lead  0.28 0.13 0.15 0.14 0.18 0.16 0.16 0.16 0.16

Reading the output. A one-standard-deviation shock to the indicator raises it by \(0.28\) at once, and it settles \(0.16\) higher for good. Sales do not respond for three periods, then jump by \(1.31\) and climb towards \(2.87\). The long-run ratio, \(2.87/0.16 \approx 18\), is the cointegrating coefficient: in the long run sales move with the indicator in the proportion fixed by \(\beta \). The three-period delay is the lag that Box and Jenkins found. Ordering sales first would change little here, because the contemporaneous correlation of the two disturbances is small; in macroeconomic VARs the ordering often matters a great deal.

Stata: var sales lead, lags(1/3), then irf create and irf graph oirf. EViews: View \(\to \) Impulse Response in a VAR object, choosing the Cholesky ordering.