Experimetrics: 9.3.1 Simulating data from a linear model

今回はExperimetrics本の9.3.1を再現してみます。先に言っておきます。この9.3.1は誤植が多いです。


9.3.1では、まず以下のような線形モデルからシミュレーションデータを生成することになってます(p.214の9.1を参照してください)。
\begin{align*}
y_{i} &= 2.0 + 0.5 x_{i} +u_{i}\;\;\;\;i=1,2,\dots,100\\
x_{i}&\sim U(0,1)\\
u_{i}&\sim N(0,1)
\end{align*}

しかし、実際には
\begin{align*}
y_{i} &= 2.0 + x_{i} +0.5 u_{i}\;\;\;\;i=1,2,\dots,100\\
x_{i}&\sim U(0,1)\\
u_{i}&\sim N(0,1)
\end{align*}
からデータが生成されたようです。

とりあえず、9.1のモデルから生成されたシミュレーションデータを使って推定を行います。使用したRコードは次の通りです。

set.seed(123)
iter<-1000
intercept<-numeric(iter)
slope<-numeric(iter)
sigma<-numeric(iter)

for (i in 1:iter) {
  # DGP (data generating process)
  n<-100
  x<-runif(n,0,1)
  u<-rnorm(n,0,1)
  y<-2+0.5*x+u
  
  fit<-lm(y~x)
  intercept[i]<-coef(fit)[1]
  slope[i]<-coef(fit)[2]
  
  res<-residuals(fit)
  sigma[i]<-sqrt(sum(res^2)/(n-2)) # or sigma[i]<-sigma(fit)
}
library(psych)
describe(intercept)
describe(slope)
describe(sigma)

シミュレーションデータをもとにOLSで推定した切片をinterceptに、傾きをslopeに、誤差項の標準誤差をsigmaに保存します。この作業を1000回繰り返すので、結果的に1000セット分の推定値を得ることになります。

結果はこちら。まずはinterceptから。1000回分の推定値の平均が2.01。データ生成に使用したinterceptが2.0ですからよく推定できてます。

   vars    n mean   sd median trimmed  mad  min max range  skew kurtosis   se
X1    1 1000 2.01 0.19   2.01    2.01 0.19 1.38 2.7  1.32 -0.02     0.14 0.01

次にslope。推定値の平均は0.49。データ生成時のslopeは0.5ですから、こちらも推定できてます。

   vars    n mean   sd median trimmed  mad   min  max range skew kurtosis   se
X1    1 1000 0.49 0.34    0.5    0.49 0.33 -0.67 1.59  2.26 0.02     0.13 0.01

最後にsigma

   vars    n mean   sd median trimmed  mad  min  max range skew kurtosis se
X1    1 1000    1 0.07      1       1 0.07 0.79 1.27  0.49 0.23     0.31  0

1000回分の推定値の分布はこんな感じです。
f:id:hiroecon:20180524083739p:plain


肝心の本の結果の再現ですが、プログラム中のDGPを次のように変更すると、結果を再現できます。

  x<-runif(n,0,1)
  u<-rnorm(n,0,1)
  y<-2+x+0.5*u