初めまして。
Qiita初投稿です。
自分は統計学や計量経済学みたいな学問に興味があり、いろいろ勉強しているのですが、勉強したことを目に見える形で残しておきたい欲が湧いてきたので、学習内容を書き留めていこうと思います。(他人が読むことを前提にしてない備忘録的な感じですが、、、)
とはいえ、まだまだわからないことだらけですし、ここに書くことも詳しい人が見ると全然正しくないこともたくさん書いちゃうと思いますが許してください、、、(いつかは完璧な記事を書きたいけどね)
今回やること
最近「実証分析のための計量経済学」(山本勲、中央経済社)を読んでて、実際に手を動かしてみたいなと思ったので、Rの「plm」というパッケージにある「Cigar」というデータセットを用いて、普通のOLS推定、LSDV推定、GLS推定の三つをやってみようと思います。Rのコードは長倉大輔さんの資料を参考にさせていただきました。
「Cigar」データセットですが、これはアメリカの46州における、1963年から1992年までのタバコの販売量とかのパネルデータです。このデータで販売量を非説明変数にとって、タバコの価格や物価でOLSをやっていきたいと思います。
本当は仮説とかきちんと立ててから分析やるべきですが、まあ今回はplmを使ってみることが目的なので、割愛します(多分、price(価格)はsales(販売量)に影響を与えているだろうなぁ)
データ読み込み
plmパッケージをインストール、読み込み。データを確認します
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をしてみましょう
resultOLS <- plm(sales~price+pop16+cpi+ndi+pimin,data = Cigar, model = "pooling")
summary(resultOLS)
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乗法は一番シンプルな形だとこんな感じになると思われます
Y_{it} = a + bX_{it} + u_{it} = a + bX_{it} +(F_i+v_{it})
パネルデータのモデルは経済主体(今回だと州)iと時点tの二つの添え字を使って表すことが一般的です。
F_iは時間によって変わらない経済主体に固有の要素を表していて、固定効果(個別効果)と呼ばれます。
これはもともと誤差u_itに含まれているものですが、この固有の要素F_iが大きいと、ある経済主体は常に大きな誤差を持ってしまうことになり、誤差の間での独立性がなくなってしまいます(共分散を持つ)
そこで、パネルデータでは誤差に含まれる州に固有の要素(固定効果)を分けて考えるのが一般的だそう。
最小2乗ダミー変数(LSDV)推定
今回のモデル式はこんな感じです
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)推定を行う。これは州ごと(観測個体)ごとのダミー変数をとって以下のようなモデル式となっています
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となるダミー変数です。このモデルだと、固定効果を調整した上で、推定結果が得られます(多分)。
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検定です。
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)を行います。(なんでかはあんまよくわからんかった)
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でこの検定はハウスマン検定と呼ばれています。
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がどういうことに向いてるやり方なのかとかふに落ちてない部分はたくさんあるので、基本から勉強したい(しなければならない)と思った。
あと、本当はパネルデータを自分で用意してやりたかったが、いいデータが見つからなかったのと、時間がかかりそうだったので今回は見送った。1個人が使えるデータは意外ときちんと探さないとないので、使えそうなデータに対してのアンテナをはっておかないとなと思ったり。
今回は計量経済学のことについて書いたが、統計学よりも、仮定を厳密に守ろうとしてる感があった。実際これからデータの分析の中でも、どんな分野を重点的に勉強していこうかは決めあぐねているところ。統計学の基本的な勉強を進めつつも、もう少し難しいことや、今回のようなデータ解析にチャレンジしていきたい。