Rでガウス＝マルコフの定理

seed <- 20250707
options(scipen = 999)

# 必要なパッケージをロード
library(ggplot2)

# ----- 1. シミュレーションの準備 -----
set.seed(seed) # 結果の再現性を確保するための乱数シード

n_sim <- 100000 # シミュレーションの繰り返し回数
n_obs <- 100 # 1回あたりのサンプルサイズ

# 真のパラメータを設定
beta0_true <- 5 # 真の切片
beta1_true <- 2 # 真の傾き
sigma_e <- 10 # 誤差項の標準偏差 (均一分散)

# 説明変数 x を生成 (シミュレーション中は固定)
x <- 1:n_obs

# 結果を保存するためのベクトルを初期化
beta1_ols <- numeric(n_sim) # OLSによる推定値を保存
beta1_alt <- numeric(n_sim) # 代替案による推定値を保存


# ----- 2. シミュレーションの実行 -----
for (i in 1:n_sim) {
  # ガウス=マルコフの仮定を満たす誤差項を生成
  # E[e]=0, Var(e)=sigma_e^2, 互いに無相関
  epsilon <- rnorm(n = n_obs, mean = 0, sd = sigma_e)

  # 観測データ y を生成
  y <- beta0_true + beta1_true * x + epsilon

  # 【推定量A: OLS (BLUE)】
  # lm()関数で最小二乗法による回帰分析を実行
  model_ols <- lm(y ~ x)
  beta1_ols[i] <- coef(model_ols)[2] # 傾きの推定値を保存

  # 【推定量B: 代替の線形不偏推定量】
  # 最初と最後のデータ点のみを使って傾きを計算
  # この推定量も線形(yの線形結合)かつ不偏(期待値がbeta1_trueになる)
  beta1_alt[i] <- (y[n_obs] - y[1]) / (x[n_obs] - x[1])
}

# ----- 3. 結果の分析 -----
cat("シミュレーション結果の要約:\n")
cat("----------------------------------\n")
cat("真の傾き (β1):", beta1_true, "\n\n")

cat("【OLS推定量】\n")
cat("  平均値:", mean(beta1_ols), " (不偏性の確認)\n")
cat("  真の価(=2)からの乖離:", abs(beta1_true - mean(beta1_ols)), "\n")
cat("  分散:  ", var(beta1_ols), " (効率性の確認)\n\n")

cat("【代替推定量】\n")
cat("  平均値:", mean(beta1_alt), " (不偏性の確認)\n")
cat("  真の価(=2)からの乖離:", abs(beta1_true - mean(beta1_alt)), "\n")
cat("  分散:  ", var(beta1_alt), " (効率性の確認)\n")
cat("----------------------------------\n")
cat(sprintf(
  "OLS推定量の分散は、代替推定量の分散の %.2f %% です。\n",
  var(beta1_ols) / var(beta1_alt) * 100
))