Chapter 8. Restricted Estimation • hansen

Load the data and create variables

We process the data with the following code:

library(hansen)
#> 
#> Attaching package: 'hansen'
#> The following object is masked from 'package:stats':
#> 
#>     influence

xx <- as.matrix(cbind(
  1,
  log(mrw1992$y60),
  log(mrw1992$invest / 100),
  log(mrw1992$pop_growth / 100 + 0.05),
  log(mrw1992$school / 100)
))

lndY <- log(mrw1992$y85) - log(mrw1992$y60)

y <- as.matrix(lndY[mrw1992$n == "Yes"])
x <- xx[mrw1992$n == "Yes", ]

Table 8.1

Using the OLS functions from Chapter 4, we can compute the slopes and standard errors with the following C++ code:

#include "00_main.h"

Mat<double> ols3_(const doubles_matrix<>& y, const doubles_matrix<>& x) {
  Mat<double> Y = as_Mat(y);
  Mat<double> X = as_Mat(x);

  int n = x.nrow();
  int k = x.ncol();

  Mat<double> beta_ols = ols_(Y, X);
  Mat<double> e_ols = Y - (X * beta_ols);
  Mat<double> Xe_ols = X.each_col() % e_ols;
  Mat<double> invX = inv(X.t() * X);
  Mat<double> V_ols = (n / (n - k)) * invX * (Xe_ols.t() * Xe_ols) * invX;
  Mat<double> se_ols = sqrt(V_ols.diag());

  Mat<double> res = join_rows(beta_ols, se_ols);
  return res;
}

[[cpp11::register]] doubles_matrix<> ols3_r_(const doubles_matrix<>& y,
                                             const doubles_matrix<>& x) {
  return as_doubles_matrix(ols3_(y, x));
}

Now we can use the function in R:

table81_1 <- ols3(y, x)

rownames(table81_1) <- c(
  "Intercept", "log(GDP 1960)", "log(1/GDP)",
  "log(n+g+delta)", "log(School)"
)
colnames(table81_1) <- c("Estimate", "Std. Error")

round(table81_1, 3)
#>                Estimate Std. Error
#> Intercept         3.022      0.718
#> log(GDP 1960)    -0.288      0.053
#> log(1/GDP)        0.524      0.105
#> log(n+g+delta)   -0.506      0.230
#> log(School)       0.231      0.065

Now we create a new function for the constrained regression in C++:

Mat<double> cns_(const doubles_matrix<>& y, const doubles_matrix<>& x,
                 const doubles_matrix<>& r) {
  Mat<double> Y = as_Mat(y);
  Mat<double> X = as_Mat(x);
  Mat<double> R = as_Mat(r);

  int n = x.nrow();
  int k = x.ncol();

  Mat<double> invX = inv(X.t() * X);
  Mat<double> iR = invX * R * inv(R.t() * invX * R) * R.t();
  Mat<double> beta_ols = ols_(Y, X);
  Mat<double> beta_cns = beta_ols - (iR * beta_ols);
  Mat<double> e_cns = Y - (X * beta_cns);
  Mat<double> Xe_cns = X.each_col() % e_cns;
  Mat<double> V_tilde = (n / (n - k + 1)) * invX * (Xe_cns.t() * Xe_cns) * invX;
  Mat<double> V_cls = V_tilde - (iR * V_tilde) - (V_tilde * iR.t()) +
    (iR * V_tilde * iR.t());
  Mat<double> se_cns = sqrt(V_cls.diag());

  Mat<double> res = join_rows(beta_cns, se_cns);
  return res;
}

[[cpp11::register]] doubles_matrix<> cns_r_(const doubles_matrix<>& y,
                                            const doubles_matrix<>& x,
                                            const doubles_matrix<>& r) {
  return as_doubles_matrix(cns_(y, x, r));
}

Now we can use the function in R:

r <- as.matrix(c(0, 0, 1, 1, 1))

table81_2 <- cns(y, x, r)

rownames(table81_2) <- rownames(table81_1)
colnames(table81_2) <- colnames(table81_1)

round(table81_2, 3)
#>                Estimate Std. Error
#> Intercept         2.457      0.430
#> log(GDP 1960)    -0.298      0.052
#> log(1/GDP)        0.501      0.090
#> log(n+g+delta)   -0.736      0.076
#> log(School)       0.235      0.064

Now we create a new function for the efficient minimum distance regression in C++:

Mat<double> emd_(const doubles_matrix<>& y, const doubles_matrix<>& x,
                 const doubles_matrix<>& r) {
  Mat<double> Y = as_Mat(y);
  Mat<double> X = as_Mat(x);
  Mat<double> R = as_Mat(r);

  int n = x.nrow();
  int k = x.ncol();

  Mat<double> invX = inv(X.t() * X);
  Mat<double> beta_ols = ols_(Y, X);
  Mat<double> e_ols = Y - (X * beta_ols);
  Mat<double> Xe_ols = X.each_col() % e_ols;
  Mat<double> V_ols = (n / (n - k)) * invX * (Xe_ols.t() * Xe_ols) * invX;

  Mat<double> V_r = V_ols * R * inv(R.t() * V_ols * R) * R.t();
  Mat<double> beta_emd = beta_ols - V_r * beta_ols;
  Mat<double> e_emd = Y - (X * beta_emd);
  Mat<double> Xe_emd = X.each_col() % e_emd;

  Mat<double> V2 = (n / (n - k + 1)) * invX * (Xe_emd.t() * Xe_emd) * invX;
  Mat<double> V_emd = V2 - V2 * R * inv(R.t() * V2 * R) * R.t() * V2;
  Mat<double> se_emd = sqrt(V_emd.diag());

  Mat<double> res = join_rows(beta_emd, se_emd);
  return res;
}

[[cpp11::register]] doubles_matrix<> emd_r_(const doubles_matrix<>& y,
                                            const doubles_matrix<>& x,
                                            const doubles_matrix<>& r) {
  return as_doubles_matrix(emd_(y, x, r));
}

Now we can use the function in R:

table81_3 <- emd(y, x, r)

rownames(table81_3) <- rownames(table81_1)
colnames(table81_3) <- colnames(table81_1)

round(table81_3, 3)
#>                Estimate Std. Error
#> Intercept         2.481      0.435
#> log(GDP 1960)    -0.297      0.052
#> log(1/GDP)        0.462      0.079
#> log(n+g+delta)   -0.711      0.071
#> log(School)       0.249      0.063

Putting all the results together:

table81 <- cbind(
  paste0(
    round(table81_1[, "Estimate"], 3), " (",
    round(table81_1[, "Std. Error"], 3), ")"
  ),
  paste0(
    round(table81_2[, "Estimate"], 3), " (",
    round(table81_2[, "Std. Error"], 3), ")"
  ),
  paste0(
    round(table81_3[, "Estimate"], 3), " (",
    round(table81_3[, "Std. Error"], 3), ")"
  )
)

rownames(table81) <- rownames(table81_1)
colnames(table81) <- c("OLS", "CNS", "EMD")

table81
#>                OLS              CNS              EMD             
#> Intercept      "3.022 (0.718)"  "2.457 (0.43)"   "2.481 (0.435)" 
#> log(GDP 1960)  "-0.288 (0.053)" "-0.298 (0.052)" "-0.297 (0.052)"
#> log(1/GDP)     "0.524 (0.105)"  "0.501 (0.09)"   "0.462 (0.079)" 
#> log(n+g+delta) "-0.506 (0.23)"  "-0.736 (0.076)" "-0.711 (0.071)"
#> log(School)    "0.231 (0.065)"  "0.235 (0.064)"  "0.249 (0.063)"