学習記録1（固定効果モデル）

初めまして。
Qiita初投稿です。

自分は統計学や計量経済学みたいな学問に興味があり、いろいろ勉強しているのですが、勉強したことを目に見える形で残しておきたい欲が湧いてきたので、学習内容を書き留めていこうと思います。（他人が読むことを前提にしてない備忘録的な感じですが、、、）

とはいえ、まだまだわからないことだらけですし、ここに書くことも詳しい人が見ると全然正しくないこともたくさん書いちゃうと思いますが許してください、、、（いつかは完璧な記事を書きたいけどね）

今回やること
---
最近[「実証分析のための計量経済学」（山本勲、中央経済社）](https://books.google.co.jp/books/about/%E5%AE%9F%E8%A8%BC%E5%88%86%E6%9E%90%E3%81%AE%E3%81%9F%E3%82%81%E3%81%AE%E8%A8%88%E9%87%8F%E7%B5%8C%E6%B8%88%E5%AD%A6.html?id=XjuMjgEACAAJ&source=kp_book_description&redir_esc=y)を読んでて、実際に手を動かしてみたいなと思ったので、Rの「plm」というパッケージにある「Cigar」というデータセットを用いて、普通のOLS推定、LSDV推定、GLS推定の三つをやってみようと思います。Rのコードは[長倉大輔さんの資料](http://user.keio.ac.jp/~nagakura/R/R_panel.pdf)を参考にさせていただきました。

「Cigar」データセットですが、これはアメリカの46州における、1963年から1992年までのタバコの販売量とかのパネルデータです。このデータで販売量を非説明変数にとって、タバコの価格や物価でOLSをやっていきたいと思います。

本当は仮説とかきちんと立ててから分析やるべきですが、まあ今回はplmを使ってみることが目的なので、割愛します（多分、price(価格)はsales(販売量)に影響を与えているだろうなぁ）

データ読み込み
---
plmパッケージをインストール、読み込み。データを確認します

```R:cigar.R
install.packages("plm")

library(plm)
data("Cigar", package = "plm")
head(Cigar, 10)
```
結果はこんな感じ

```
   state year price    pop  pop16   cpi       ndi sales pimin
1      1   63  28.6 3383.0 2236.5  30.6  1558.305  93.9  26.1
2      1   64  29.8 3431.0 2276.7  31.0  1684.073  95.4  27.5
3      1   65  29.8 3486.0 2327.5  31.5  1809.842  98.5  28.9
4      1   66  31.5 3524.0 2369.7  32.4  1915.160  96.4  29.5
5      1   67  31.6 3533.0 2393.7  33.4  2023.546  95.5  29.6
6      1   68  35.6 3522.0 2405.2  34.8  2202.486  88.4  32.0

```


OLS推定
---

右から第二列目のsales（一人当たりのタバコ販売量(箱)）を非説明変数として使います。まずは普通にOLSをしてみましょう

```R:cigar.R
resultOLS <- plm(sales~price+pop16+cpi+ndi+pimin,data = Cigar, model = "pooling")
summary(resultOLS)
```

```output
Coefficients:
               Estimate  Std. Error  t-value  Pr(>|t|)    
(Intercept)  1.3770e+02  2.0047e+00  68.6903 < 2.2e-16 ***
price       -1.5305e+00  1.1527e-01 -13.2770 < 2.2e-16 ***
pop16       -9.2129e-04  2.1193e-04  -4.3472 1.481e-05 ***
cpi          1.5581e-01  8.0208e-02   1.9426   0.05226 .  
ndi          6.1766e-03  5.7815e-04  10.6834 < 2.2e-16 ***
pimin        5.8099e-01  1.2653e-01   4.5919 4.794e-06 ***
```
ふむふむ、price(価格)、pop16(16才以上の人口)、ndi(一人当たり可処分所得)、pimin(隣の州のタバコ一箱の最低価格)のp値が非常に小さいですね。cpi(アメリカの消費者物価指数)だけがp値が大きくなっているという感じです。

パネルデータ分析の注意点
---
ここで、突然ですが最小2乗法は一番シンプルな形だとこんな感じになると思われます

```math
Y_{it} = a + bX_{it} + u_{it} = a + bX_{it} +(F_i+v_{it})
```
パネルデータのモデルは経済主体（今回だと州）iと時点tの二つの添え字を使って表すことが一般的です。

F_iは時間によって変わらない経済主体に固有の要素を表していて、固定効果（個別効果）と呼ばれます。

これはもともと誤差u_itに含まれているものですが、この固有の要素F_iが大きいと、ある経済主体は常に大きな誤差を持ってしまうことになり、誤差の間での独立性がなくなってしまいます（共分散を持つ）

そこで、パネルデータでは誤差に含まれる州に固有の要素（固定効果）を分けて考えるのが一般的だそう。

最小2乗ダミー変数(LSDV)推定
---
今回のモデル式はこんな感じです

```math
sales_{it} = \alpha_i + \beta_1 price_{it} + \beta_2 pop16_{it} + ... + \beta_5pimin_{it} + \varepsilon_{it} \\
i = 1,2,...46\quad t =1,2,...30
```

iは州の数、tは年(1963年から1992年)を表します。固定効果に当たるのは今回だと、α_iですね。

普通のOLSでは、α_1=α_2=...=α_Nとして推定を行っていますが、この仮定はどの州もタバコの販売量に与える効果は同じであるということ（つまりその州に固有の効果はない、固定効果はないこと）を意味します。実際はタバコの販売量に与える影響は州ごとに違うはずなので、この仮定はよくないですね。そこでα_1=α_2=...=α_Nという仮定をおかない、最小2乗ダミー変数(LSDV)推定を行う。これは州ごと（観測個体）ごとのダミー変数をとって以下のようなモデル式となっています

```math
sales_{it} = \alpha_1D1_{it}+\alpha_2D2_{it}+...+\alpha_{46}D46_{it}+\quad\beta_1 price_{it} + \beta_2 pop16_{it} + ... + \beta_5pimin_{it} + \varepsilon_{it}
```

D1_iはi=1の時D1_i=1となり、iが1以外の時はD1_i=0となるダミー変数です。このモデルだと、固定効果を調整した上で、推定結果が得られます（多分）。

```R:cigar.R
resultLSDV <- plm(sales~price+pop16+cpi+ndi+pimin,data = Cigar,model = "within")
summary(resultLSDV)
```
結果はこちら

```
Coefficients:
         Estimate  Std. Error  t-value Pr(>|t|)    
price -0.65718025  0.07550568  -8.7037  < 2e-16 ***
pop16  0.00131517  0.00057258   2.2969  0.02178 *  
cpi    0.95780088  0.05135959  18.6489  < 2e-16 ***
ndi   -0.00494601  0.00045873 -10.7819  < 2e-16 ***
pimin  0.17858026  0.08558861   2.0865  0.03712 * 
```

pop16やpiminのp値は大きくなり、cpiのp値は小さくなりました。

ここで本当に固定効果はないα_1=α_2=...=α_Nという仮定は正しくないのか検定を行ってみます
帰無仮説：α_1=α_2=...=α_NのF検定です。

```R:cigar.R
pFtest(resultLSDV,resultOLS)
```
```
	F test for individual effects

data:  sales ~ price + pop16 + cpi + ndi + pimin
F = 97.74, df1 = 45, df2 = 1329, p-value < 2.2e-16
alternative hypothesis: significant effects
```
p値をみると帰無仮説は棄却されます。よって、「固定効果はある」とわかります。

最後に一般化最小2乗法(GLS)を行います。（なんでかはあんまよくわからんかった）

```R:cigar.R
resultGLS <- plm(sales~price+pop16+cpi+ndi+pimin,data = Cigar,model = "random")
summary(resultGLS)
```

```
Coefficients:
               Estimate  Std. Error z-value Pr(>|z|)    
(Intercept)  1.2175e+02  3.5213e+00 34.5766  < 2e-16 ***
price       -6.7400e-01  7.6010e-02 -8.8672  < 2e-16 ***
pop16        8.5822e-04  4.8662e-04  1.7637  0.07779 .  
cpi          9.3235e-01  5.1581e-02 18.0752  < 2e-16 ***
ndi         -4.5519e-03  4.5921e-04 -9.9125  < 2e-16 ***
pimin        1.7880e-01  8.6033e-02  2.0783  0.03769 *  
```

推定結果は得られましたが、GLSは固定効果と説明変数に相関があると、得られた推定量の一致性（サンプルサイズを大きくすれば、推定量がある値に収束する性質）がなくなってしまいます。

これはあまりいい推定量とはいえないので、固定効果と説明変数に相関があるかどうかを検定します。

帰無仮説：相関＝0でこの検定はハウスマン検定と呼ばれています。

```R:cigar.R
phtest(resultLSDV, resultGLS)
```

```
	Hausman Test

data:  sales ~ price + pop16 + cpi + ndi + pimin
chisq = 182.95, df = 5, p-value < 2.2e-16
alternative hypothesis: one model is inconsistent
```
p値より、帰無仮説は棄却されるので、相関は0ではなく、固定効果と説明変数に相関があります。よって、GLSを使うための仮定が満たされていないので、GLSの結果ではなく、LSDVの結果をみることにします。

感想
---
いろいろごまかしながら記事を書いたし、数学的な部分が一切わかっていないのはなんとかしたいと思う。
あとは固定効果が誤差間の独立性をなくす部分とか、GLSがどういうことに向いてるやり方なのかとかふに落ちてない部分はたくさんあるので、基本から勉強したい（しなければならない）と思った。

あと、本当はパネルデータを自分で用意してやりたかったが、いいデータが見つからなかったのと、時間がかかりそうだったので今回は見送った。１個人が使えるデータは意外ときちんと探さないとないので、使えそうなデータに対してのアンテナをはっておかないとなと思ったり。

今回は計量経済学のことについて書いたが、統計学よりも、仮定を厳密に守ろうとしてる感があった。実際これからデータの分析の中でも、どんな分野を重点的に勉強していこうかは決めあぐねているところ。統計学の基本的な勉強を進めつつも、もう少し難しいことや、今回のようなデータ解析にチャレンジしていきたい。