R语言带系数约束的回归及约束性2SLS回归实现求助

Great question! You’re right that rolling your own constrained OLS from scratch isn’t the most efficient approach—R has dedicated tools for both constrained least squares and constrained 2SLS that will save you time and reduce errors. Let’s walk through both scenarios step by step.

Your problem is a classic constrained least squares problem: minimize the sum of squared residuals while enforcing two linear constraints:

1. β₁ + β₂ + β₃ = 1
2. β₁ ≥ 0, β₂ ≥ 0, β₃ ≥ 0

This is a quadratic programming problem (since the sum of squares is a quadratic function), and the quadprog package is perfect for this—it’s optimized for exactly these kinds of constrained optimization tasks.

Step-by-Step Implementation #

First, install and load the package:

install.packages("quadprog")
library(quadprog)

Let’s simulate some sample data to test with:

set.seed(123) # For reproducibility
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
Y <- 0.3*x1 + 0.5*x2 + 0.2*x3 + rnorm(n, 0, 0.1) # True coefficients meet your constraints

Now construct the matrices needed for quadratic programming:

• Dmat: The cross-product matrix of your predictors (X'X), which defines the quadratic part of the objective function
• dvec: The cross-product of predictors and the response (X'Y), which defines the linear part
• Amat & bvec: These define your constraints. We’ll encode the sum-to-1 constraint as two inequality constraints (β₁+β₂+β₃ ≥ 1 and -(β₁+β₂+β₃) ≥ -1) to turn it into an equality, plus non-negativity constraints for each coefficient.

X <- cbind(x1, x2, x3)
Dmat <- t(X) %*% X
dvec <- t(Y) %*% X

# Build constraint matrix:
# Rows 1-2: Enforce β₁+β₂+β₃ = 1
# Rows 3-5: Enforce β₁, β₂, β₃ ≥ 0
Amat <- cbind(rep(1, 3), rep(-1, 3), diag(3))
bvec <- c(1, -1, 0, 0, 0)

# Solve the quadratic program: meq=2 means first 2 constraints are equalities
solution <- solve.QP(Dmat, dvec, Amat, bvec, meq = 2)
constrained_beta <- solution$solution
names(constrained_beta) <- c("x1", "x2", "x3")

# Check results
constrained_beta
sum(constrained_beta) # Should equal 1
all(constrained_beta >= 0) # Should be TRUE

For constrained 2SLS, we need to combine the two-stage least squares logic with the same coefficient constraints. There are two reliable approaches here:

Approach 1: First-Stage Fitting + Constrained OLS #

This is the more intuitive method:

1. Run first-stage regressions for any endogenous predictors using your instrumental variables
2. Replace endogenous predictors with their fitted values
3. Run constrained OLS on the transformed predictors (using the same quadprog method as above)

Example Implementation #

Let’s assume x2 is endogenous, with instrumental variables z1 and z2:

# Simulate endogenous variable and instruments
set.seed(123)
z1 <- rnorm(n)
z2 <- rnorm(n)
x2 <- 0.4*z1 + 0.3*z2 + rnorm(n, 0, 0.1) # x2 correlates with instruments
Y <- 0.3*x1 + 0.5*x2 + 0.2*x3 + rnorm(n, 0, 0.1) # Y with endogenous x2

# Step 1: First-stage regression for endogenous x2
first_stage <- lm(x2 ~ x1 + x3 + z1 + z2)
x2_hat <- predict(first_stage) # Fitted values of x2

# Step 2: Constrained OLS on Y ~ x1 + x2_hat + x3
X_2sls <- cbind(x1, x2_hat, x3)
Dmat_2sls <- t(X_2sls) %*% X_2sls
dvec_2sls <- t(Y) %*% X_2sls

# Reuse the same constraint setup as before
Amat_2sls <- cbind(rep(1, 3), rep(-1, 3), diag(3))
bvec_2sls <- c(1, -1, 0, 0, 0)

solution_2sls <- solve.QP(Dmat_2sls, dvec_2sls, Amat_2sls, bvec_2sls, meq = 2)
constrained_2sls_beta <- solution_2sls$solution
names(constrained_2sls_beta) <- c("x1", "x2", "x3")

constrained_2sls_beta

Approach 2: Direct Constrained Optimization of 2SLS Objective #

If you want to avoid separate first-stage fitting, you can directly optimize the 2SLS objective function (minimizing the projected residual sum of squares) with your constraints using constrOptim:

# Define the 2SLS objective function
objective_2sls <- function(beta, Y, X, Z) {
  P_Z <- Z %*% solve(t(Z) %*% Z) %*% t(Z) # Projection matrix for instruments
  residuals <- Y - X %*% beta
  sum(t(residuals) %*% P_Z %*% residuals) # Projected sum of squares
}

# Instrument matrix Z includes exogenous predictors + instruments
Z <- cbind(x1, x3, z1, z2)
X <- cbind(x1, x2, x3)

# Initial guess: Unconstrained 2SLS coefficients
unconstrained_2sls <- lm(Y ~ x1 + x2 + x3 - 1, instrument = ~ z1 + z2 + x1 + x3)
init_beta <- coef(unconstrained_2sls)

# Define constraints: β₁+β₂+β₃=1 and all β≥0
ui <- rbind(c(1, 1, 1), c(-1, -1, -1), diag(3)) # Constraint matrix
ci <- c(1, -1, 0, 0, 0) # Lower bounds for constraints

# Run constrained optimization
result_2sls <- constrOptim(
  theta = init_beta,
  f = objective_2sls,
  grad = NULL, # Let R compute numerical gradients
  ui = ui,
  ci = ci,
  Y = Y,
  X = X,
  Z = Z
)

# Get constrained coefficients
result_2sls$par

Key Notes #

• The quadprog approach is faster for constrained OLS since it’s specialized for quadratic problems.
• For constrained 2SLS, both methods work, but Approach 1 is easier to debug, while Approach 2 is more direct.

内容的提问来源于stack exchange，提问作者Brad G.