すごく基本的な線形回帰のモンテカルロ・シミュレーション
知っている人にとってはすごく基本的なことだけれど、復習がてら。
ステップ
コード
まず、Aという「入れ物」を作る。
A= rep(NA,500)
乱数Xを用意する。
x = rnorm(25,mean=2,sd=1)
500回ループを回す。
一行目は、条件付き変数YをXから作成。
二行目で、XをYに回帰させる推定。
3行目でbetaをAに保存する。
for(i in 1:500){ y = rnorm(25, mean=(3*x+2), sd=1) beta = lm(y~x) A[i] = beta$coef[2] }
推定されたの平均と分散を計算。
Abar = mean(A) varA = var(A)
結果は、平均が3.026461, 分散が0.0710012。
結果のヒストグラムは以下。真の値である3を中心に、散っている感じ。
結果
線形回帰の仮定が満たされているデータ生成仮定(DGP)のシミュレーションをしたわけですが、
この結果は果たして妥当なのか?つまり、真の値が満たされているのか?
せっかくなので、復習がてら手動でt検定。
帰無仮説は
T統計量は
n = length(x) # サンプル数 k = length(beta$coefficients) # パラメータの数 seA = sd(A)/sqrt(n) #標準誤差の計算 tstat = (Abar-3)/seA # t統計量の計算
t統計量は、自由度のt分布に従うので、その分布で検定。
このp値が、両側検定で5%なら真の値が3であるという仮説は棄却される。つまり結果が0.95より大きいか、0.05より小さければ棄却。
pt(tstat,n-k)
結果は0.6878834。なので、帰無仮説は棄却されない。
回帰分析の仮定が満たされてるデータを生成して回帰分析(最小二乗法)した結果、バイアスのない値が推定されることがわかりました。
すごく基本だけど、久しぶりにやり直すと勉強になる。
元ネタは、Econometrics in Rという本のPDFの付録から。文字通り「Rで計量経済学」ですね。
https://ocw.mit.edu/courses/economics/14-381-statistical-method-in-economics-fall-2013/study-materials/MIT14_381F13_EcnomtrisInR.pdf