Rで最小二乗推定量の不偏性と一致性

# --- 0. 準備: 必要なライブラリの読み込み ---
library(ggplot2)
library(dplyr)
library(tidyr)

result_output <- "D:/result_output"

# --- 1. シミュレーションの基本設定 ---
seed <- 20250704
n_sim <- 2000 # 各シナリオ・サンプルサイズでのシミュレーション回数
true_beta0 <- 2.0 # 真の切片 (α)
true_beta1 <- 0.5 # 真の傾き (β)

# テストするサンプルサイズ
sample_sizes <- c(50, 200, 1000, 4000)

# 系列相関のパラメータ
rho <- 0.8 # AR(1)モデルの自己相関係数 u_t = rho * u_{t-1} + e_t

# 結果を格納するための空のリスト
results_list <- list()

# --- 2. シミュレーションの実行 ---

set.seed(seed)

# サンプルサイズごとにループ
for (n in sample_sizes) {
  # シミュレーション回数だけループ
  for (i in 1:n_sim) {
    # --- パターン1: 誤差項が正規分布でない ---
    x1 <- rnorm(n, mean = 10, sd = 3)
    u1 <- rt(n, df = 3)
    y1 <- true_beta0 + true_beta1 * x1 + u1
    fit1 <- lm(y1 ~ x1)

    # --- パターン2: 系列相関あり & 説明変数が厳密に外生的 ---
    x2 <- rnorm(n, mean = 10, sd = 3)
    u2 <- arima.sim(model = list(ar = rho), n = n, sd = 1)
    y2 <- true_beta0 + true_beta1 * x2 + u2
    fit2 <- lm(y2 ~ x2)

    # --- パターン3: 系列相関あり & 説明変数にラグ変数を含む ---
    u3 <- numeric(n + 1)
    y3 <- numeric(n + 1)
    e3 <- rnorm(n + 1, sd = 1)
    y3[1] <- true_beta0 / (1 - true_beta1)
    u3[1] <- e3[1] / sqrt(1 - rho^2)
    for (t in 2:(n + 1)) {
      u3[t] <- rho * u3[t - 1] + e3[t]
      y3[t] <- true_beta0 + true_beta1 * y3[t - 1] + u3[t]
    }
    y_t <- y3[2:(n + 1)]
    y_lag <- y3[1:n]
    fit3 <- lm(y_t ~ y_lag)

    # 結果をリストに追加
    results_list[[length(results_list) + 1]] <- data.frame(
      sample_size = n,
      simulation_id = i,
      scenario = c(
        "1. 誤差項が非正規分布",
        "2. 系列相関あり (説明変数が外生的)",
        "3. 系列相関あり (ラグ変数を含む)"
      ),
      estimated_beta1 = c(coef(fit1)[2], coef(fit2)[2], coef(fit3)[2])
    )
  }
  # 進捗表示
  # cat("サンプルサイズ n =", n, "のシミュレーションが完了しました。\n")
}

# リストを一つのデータフレームに変換
results_df <- bind_rows(results_list)

# シミュレーション結果の保存
setwd(result_output)
saveRDS(object = results_df, file = "ols_unbiasedness_consistency.rds")

# tbl オブジェクトの読み込み
setwd(result_output)
results_df <- readRDS("ols_unbiasedness_consistency.rds")

# --- 3. 結果の確認（不偏性と一致性） ---

# 3-1. 不偏性の確認 (推定量の平均値)
unbiased_check <- results_df %>%
  group_by(scenario, sample_size) %>%
  summarise(mean_beta1_hat = mean(estimated_beta1), .groups = "drop") %>%
  mutate(bias = mean_beta1_hat - true_beta1)

cat("--- 不偏性の確認 (真のβ1 = ", true_beta1, ") ---\n")
unbiased_check %>% knitr::kable()