在R中生成随机样本以对比GLS与OLS估计量渐近性的问题

需求背景 #

需要对比GLS（广义最小二乘）与OLS（普通最小二乘）估计量的渐近性，计划生成6000个随机样本：

• 样本量：50、100两种
• $\gamma_0$参数取值：0、0.1、0.5三种
• 上述6种组合各重复1000次

所用模型：
$$y_i=\alpha+\betax_i+u_i$$
其中 $u_i=e_ih(x_i)^{1/2}$，$h(x)=exp(\gamma_0*x)$

原代码问题分析 #

1. 参数引用错误：循环中用gamma0[k]调用参数，但k是$\gamma_0$的实际取值（如0、0.1），不是索引，导致h计算出错为空
2. 样本存储覆盖：用u[[i]]存储结果会覆盖同一样本量下不同$\gamma_0$参数的样本
3. 缺少重复循环：未实现1000次重复逻辑，仅生成了单次样本
4. 未生成因变量y：后续估计OLS/GLS需要用到y，但原代码未定义

修复后的完整代码 #

# 模拟设置
N <- 1000  # 每个组合的重复次数
sample_sizes <- c(50, 100)
gamma_vals <- c(0, 0.1, 0.5)
true_params <- list(alpha = 1, beta = 1)

# 初始化结果存储矩阵：N行（重复次数）×6列（6种组合）
alpha_ols <- matrix(NA, nrow = N, ncol = 6)
beta_ols <- matrix(NA, nrow = N, ncol = 6)
alpha_gls <- matrix(NA, nrow = N, ncol = 6)
beta_gls <- matrix(NA, nrow = N, ncol = 6)

# 组合索引：对应6列的顺序
comb_idx <- 1

# 三层循环：样本量 → gamma参数 → 重复次数
for (n in sample_sizes) {
  for (gamma in gamma_vals) {
    for (rep in 1:N) {
      # 生成自变量x和随机误差e
      x <- rnorm(n, mean = 0, sd = 1)
      e <- rnorm(n, mean = 0, sd = 1)
      
      # 计算异方差项h和扰动项u
      h <- exp(gamma * x)
      u <- e * sqrt(h)
      
      # 生成因变量y
      y <- true_params$alpha + true_params$beta * x + u
      
      # 计算OLS估计
      ols_model <- lm(y ~ x)
      alpha_ols[rep, comb_idx] <- coef(ols_model)[1]
      beta_ols[rep, comb_idx] <- coef(ols_model)[2]
      
      # 计算GLS估计（加权最小二乘，权重为1/h）
      gls_model <- lm(y ~ x, weights = 1/h)
      alpha_gls[rep, comb_idx] <- coef(gls_model)[1]
      beta_gls[rep, comb_idx] <- coef(gls_model)[2]
    }
    # 切换到下一个组合的列
    comb_idx <- comb_idx + 1
  }
}

代码说明 #

• 组合索引comb_idx：按(50,0)→(50,0.1)→(50,0.5)→(100,0)→(100,0.1)→(100,0.5)的顺序，将每个组合的1000次估计结果存入对应列
• GLS实现：通过加权最小二乘lm(..., weights = 1/h)近似GLS，因为已知异方差结构$Var(u_i)=h(x_i)$
• 结果存储：所有估计值按重复次数和组合分别存入对应矩阵，后续可直接计算均值、方差等统计量分析渐近性

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