［R］異常値を含むダミー変数を用いた回帰直線の推定（「44の例題で学ぶ計量経済学」（オーム社）pp.226-227）

> dtf <- read.csv("table5_1.csv", header = TRUE)
> n <- nrow(dtf)
> xxi <- dtf$xxi
> ddi <- dtf$ddi
> yyi <- dtf$yyi
> k <- 3
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- matrix(c(rep(1.0, n), xxi, ddi), ncol = 3)
> mxbh <- solve(t(mxxx) %*% mxxx) %*% t(mxxx) %*% mxy
> mxyh <- mxxx %*% mxbh
> mxu <- mxy - mxyh
> ssy2 <- sum(t(mxy) %*% mxy) - sum(yyi) ^ 2 / n
> ssyh2 <- sum(t(mxyh) %*% mxyh) - sum(yyi) ^ 2 / n
> sse2 <- as.vector(t(mxu) %*% mxu)
> s2 <- sse2 / (n - k)
> s <- sqrt(s2)
> mxvv <- s2 * solve(t(mxxx) %*% mxxx)
> mxs <- sqrt(diag(mxvv))
> rr2 <- ssyh2 / ssy2
> tval <- as.vector(mxbh) / as.vector(mxs)
> # 推定された回帰モデルのパラメーター α^,β^,γ^
> print(as.vector(mxbh))
[1]  0.5  1.1 15.8
> # 決定係数 R^2
> print(rr2)
[1] 0.9936886
> # 各パラメーターの t値
> print(tval)
[1]  0.4303315  5.1854497 13.8309443
> # 図の出力
> png("fig5_1.png", width = 512, height = 512)
> plot(xxi, yyi, asp = 1., type = "n", xlim = c(0, 12), ylim = c(0, 30), xlab = "X", ylab = "Y")
> x <- seq(0, 12, length = 100)
> lines(x, mxbh[1] + mxbh[2] * x, lty = "solid")
> lines(x, mxbh[1] + mxbh[2] * x + mxbh[3], lty = "dashed")
> points(xxi, yyi, pch = 20)
> dev.off()
null device 
          1