データ分析メモと北欧生活

旧Untitled Note. データ分析、計量経済・統計とR、水産管理、英語勉強、海外生活などについて備忘録や自分の勉強のOutputの場所として

MENU

すごく基本的な線形回帰のモンテカルロ・シミュレーション

簡単なモンテカルロ法を、R言語でやります。

知っている人にとってはすごく基本的なことだけれど、復習がてら。

ステップ

  1. Xという乱数を作る。
  2. Xに依存した乱数Yを作成する。
  3. 回帰分析をして、 \betaを推定する。
  4. 以上のプロセスを500回繰り返す。
    • 500個の推定された \betaから、分布を作る。

コード

まず、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]
      }

推定された \betaの平均と分散を計算。

Abar = mean(A)
varA = var(A)

結果は、平均が3.026461, 分散が0.0710012。

結果のヒストグラムは以下。真の値である3を中心に、散っている感じ。

f:id:keita43a:20180502033545p:plain

結果

線形回帰の仮定が満たされているデータ生成仮定(DGP)のシミュレーションをしたわけですが、
この結果は果たして妥当なのか?つまり、真の値 \beta = 3が満たされているのか?

せっかくなので、復習がてら手動でt検定。
帰無仮説 H_0: \beta = 3

T統計量は

 t = \frac{\hat{\beta} - \beta_0}{se(\hat{\beta})}

n = length(x) # サンプル数
k = length(beta$coefficients) # パラメータの数

seA = sd(A)/sqrt(n) #標準誤差の計算

tstat = (Abar-3)/seA # t統計量の計算

t統計量は、自由度 n-kの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