R语言OLS矩阵形式蒙特卡洛模拟r×2参数矩阵存储问题咨询

• 自变量矩阵未添加截距列：矩阵形式OLS估计带截距的系数时，需要给x补充一列全1的截距项，否则运算仅会返回斜率估计值，无法得到β₀结果
• 存储容器维度不匹配：原代码初始化的ols是一维向量，仅能存储单值，需要调整为nreps × 2的矩阵，分别存每次模拟的β₀、β₁估计值
• 参数赋值不匹配：你描述的真实参数为β₀=1、β₁=-1，原代码中参数赋值相反，修正代码已同步调整

set.seed(123)
n = 20
nreps = 10000
# 按你描述的真实参数赋值
beta_0 = 1 
beta_1 = -1 

# 初始化nreps行2列的矩阵存储结果，两列分别对应β0、β1估计值
ols = matrix(nrow = nreps, ncol = 2)
colnames(ols) = c("beta0_hat", "beta1_hat") # 可选，给列加标注方便后续调用

## 蒙特卡洛模拟循环
for (r in 1:nreps) {
  u = rnorm(n, mean = 0, sd = 1) 
  x = rnorm(n, mean = 1, sd = 2^2) 
  # 给x添加截距项列（全1列）
  x_mat = cbind(1, x)
  y = beta_0 + beta_1*x + u

  beta_hat = solve(t(x_mat) %*% x_mat) %*% (t(x_mat) %*% y)
  # 将本次估计的两个系数存入ols的第r行
  ols[r, ] = beta_hat
}

# 统计结果示例
ols_average = colMeans(ols)
ols_sd = apply(ols, 2, sd)
# 分别绘制两个系数的密度图
par(mfrow = c(1,2))
plot(density(ols[,1]), main = "β0估计值密度")
plot(density(ols[,2]), main = "β1估计值密度")

运行代码后ols就是你需要的nreps × 2维结果矩阵，每一行对应一次蒙特卡洛模拟的两个系数估计值，你可以直接基于该矩阵做后续的蒙特卡洛性质检验。

内容的提问来源于stack exchange，提问作者arhe