スライド

平成 26 年度計量経済分析演習
「 R による計量経済分析 (3)」
.
∼ R による重回帰分析 ∼
原尚幸
.
新潟大・経済
http://www.econ.niigata-u.ac.jp/˜hara/grad-ecm/
[email protected]
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
1 / 56
Contents
1
都道府県別コンビニ数のデータ
記述統計
回帰分析
2
物価と為替レートの関係
3
賃金関数の推定
4
たばこの価格と消費量のデータ
5
レポート課題
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
2 / 56
Outline
1
都道府県別コンビニ数のデータ
記述統計
回帰分析
2
物価と為替レートの関係
3
賃金関数の推定
4
たばこの価格と消費量のデータ
5
レポート課題
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
3 / 56
Outline
1
都道府県別コンビニ数のデータ
記述統計
回帰分析
2
物価と為替レートの関係
3
賃金関数の推定
4
たばこの価格と消費量のデータ
5
レポート課題
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
4 / 56
データのダウンロード
講義のサイトからデータ”cv.csv” をダウンロード
して Rdata に保存
http://www.econ.niigata-u.ac.jp/˜hara/grad-ecm/
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
5 / 56
データと分析の目的
都道府県別のコンビニの数とその要因と
考えられる変数のデータ (2005 年)
このデータを用いてコンビニ件数を説明
する最適なモデルの探索
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
6 / 56
データに含まれる変数
car : 1000 人あたりの自家用車の台数
conv : コンビニ件数
divorce : 離婚率
gpp : 県内総生産額 (億円)
gsm : 総合スーパーの件数 (イオンなど )
income : 1ヶ月の平均月収 (勤労者世帯のみ)
nm_m_30 : 未婚者割合 (30∼34 歳男性)
nm_f_30 : 未婚者割合 (30∼34 歳女性)
pop : 人口 (1000 人)
sm : 専門スーパーの件数 (ウオロク , ニトリなど )
w_pop : 生産年齢人口割合 (15∼64 歳)
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
7 / 56
データの読み込み
データをオブジェクト cvdata に読み込む
> cvdata <- read.csv("D:/Rdata/cv.csv")
pop, gpp, conv のデータを cvdata3 に代入
> cvdata3 <- cvdata[,c("pop","gpp","conv")]
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
8 / 56
基礎統計量の計算
gpp の最小値, 最大値, 平均値, 中央値, 四分位点
を計算
> gpp <- cvdata$gpp
最小値：min(gpp)
最大値：max(gpp)
平均値：mean(gpp)
中央値：median(gpp)
四分位点：quantile(gpp)
これらをまとめて計算
> summary(gpp)
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
9 / 56
基礎統計量の計算
gpp の分散・標準偏差を計算
不偏分散：var(gpp)
標準偏差：sd(gpp)
データを X1 , . . . , Xn としたときに,
不偏分散 var, 標準偏差 sd はそれぞれ
1 ∑
¯ 2,
(Xi − X)
n − 1 i=1
n
H. Hara (Niigata U.)
R による計量経済分析
1 ∑
¯ 2
(Xi − X)
n − 1 i=1
n
Nov 11-, 2014
10 / 56
基礎統計量の計算
pop, gpp, conv の基本統計量を計算
> summary(cvdata3)
pop
gpp
Min.
: 607
Min.
: 20057
1st Qu.: 1164
1st Qu.: 37327
Median : 1753
Median : 59248
Mean
: 2718
Mean
:109823
3rd Qu.: 2762
3rd Qu.:104927
Max.
:12577
Max.
:922694
H. Hara (Niigata U.)
R による計量経済分析
conv
Min.
: 142.0
1st Qu.: 366.5
Median : 493.0
Mean
: 909.3
3rd Qu.: 896.5
Max.
:5453.0
Nov 11-, 2014
11 / 56
基本統計量の計算
pop, gpp, conv の分散と共分散を計算
> var(cvdata3)
pop
gpp
conv
pop
6739226
356392802
2518575.1
gpp 356392802 21998634289 141292516.9
conv
2518575
141292517
991169.4
分散共分散行列
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
12 / 56
基本統計量の計算
pop, gpp の相関係数を計算
> cor(pop,gpp)
[1] 0.9256057
pop, gpp, conv の相関係数をまとめて計算
> cor(cvdata3)
pop
gpp
conv
pop 1.0000000 0.9256057 0.9744868
gpp 0.9256057 1.0000000 0.9568577
conv 0.9744868 0.9568577 1.0000000
この行列を相関行列という.
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
13 / 56
Outline
1
都道府県別コンビニ数のデータ
記述統計
回帰分析
2
物価と為替レートの関係
3
賃金関数の推定
4
たばこの価格と消費量のデータ
5
レポート課題
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
14 / 56
単回帰モデル
コンビニ数を被説明変数, 人口を説明変数とした
単回帰モデルを OLS で推定する.
convi = α + βpopi + ui
コンビニ数は人口に比例することが予想できる.
⇒β>0
まずは縦軸：コンビニ数, 横軸：人口でプロット .
> plot(conv ˜ pop, data = cvdata)
>
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
15 / 56
3000
0
1000
2000
conv
4000
5000
推定結果
2000
4000
6000
8000
10000
12000
pop
単回帰モデルがよくあてはまりそうだ
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
16 / 56
単回帰分析
OLS で推定
> sr <- lm(conv ˜ pop, data=cvdata)
> result.sr <- summary(sr)
> result.sr
推定結果
> result.sr$coefficients
Estimate Std. Error
t value
Pr(>|t|)
(Intercept) -106.6232310 47.98637570 -2.221948 3.135850e-02
pop
0.3737187 0.01283136 29.125422 7.649012e-31
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
17 / 56
推定結果
散布図に推定直線を引く
> abline(sr)
R2
3000
0
1000
2000
conv
4000
5000
> result.sr$r.squared
[1] 0.9496244
2000
4000
6000
8000
10000
12000
pop
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
18 / 56
推定結果
残差二乗和の計算
> rd <- sr$residuals
> rss <- sum(rdˆ2)
> rss
[1] 2296813
>
rd : 残差ベクトル (残差の系列)
sum(rdˆ2) : 残差二乗和
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
19 / 56
重回帰モデル
コンビニ数を被説明変数とし
人口
自家用自動車台数
女性未婚者割合
を説明変数とした重回帰モデルを OLS で推定.
convi = α + β2 popi + β3 cari + β4 nm_f_30i + ui
β2 > 0, β3 > 0, β4 > 0
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
20 / 56
重回帰分析
OLS で推定
> mr <- lm(conv ˜ pop + car + nm_f_30, data=cvdata)
> result.mr <- summary(mr)
> result.mr
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
21 / 56
推定結果
推定結果の見方は単回帰モデルのときと同じ
回帰係数
> result.mr$coefficients
Estimate Std. Error
t value
Pr(>|t|)
(Intercept) -2355.8703190 652.1889897 -3.612251 7.890445e-04
pop
0.3814832
0.0161185 23.667415 2.744634e-26
car
2.1778204
0.7345106 2.964995 4.922102e-03
nm_f_30
47.1853806 14.0482351 3.358812 1.648456e-03
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
22 / 56
推定結果
result.mr に含まれる情報
> attributes(result.mr)
$names
[1] "call"
"terms"
[4] "coefficients" "aliased"
[7] "df"
"r.squared"
[10] "fstatistic"
"cov.unscaled"
"residuals"
"sigma"
"adj.r.squared"
$class
[1] "summary.lm"
fstatistic : H0 : β1 = β2 = β3 = 0 の検定の F 統計量
adj.r.squared : 自由度調整済 R2
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
23 / 56
演習 1：線形制約の検定
演習 1
重回帰モデル
convi = α + β2 popi + β3 cari + β4 nm_f_30i + ui
において,
H0 : β3 = β4 = 0
の線形制約の検定を行え.
F 統計量 F の F (a, b) 分布における p 値の計算は
1 − pf(F, a, b)
.
pf() : 下側パーセント点の計算
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
24 / 56
パッケージ ”AER”のインストール
Kleiber and Zeileis の教科書の付属のパッケージ
パッケージをインストールして読み込む
> install.packages("AER")
> library("AER")
R には AER に限らずさまざまなパッケージが用意
されている.
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
25 / 56
線形制約の検定
線形制約の検定は AER に組み込まれている
linear.hypothesis()
という関数を用いてもできる.
> lh <- linear.hypothesis(mr,c("car=0","nm_f_30=0"))
linear.hypothesis(H1 の推定結果, 線形制約)
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
26 / 56
線形制約の検定
> lh
Linear hypothesis test
Hypothesis:
car = 0
nm_f_30 = 0
Model 1: restricted model
Model 2: conv ˜ pop + car + nm_f_30
Res.Df
RSS Df Sum of Sq
F
Pr(>F)
1
45 2296813
2
43 1783227 2
513586 6.1922 0.004333 **
--Signif. codes: 0 ‘ *** ’ 0.001 ‘ ** ’ 0.01 ‘ * ’ 0.05 ‘ . ’ 0.1 ‘
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
’1
27 / 56
AIC, BIC の計算
AIC, BIC は回帰分析の結果からそれぞれ
AIC(mr), BIC(mr) によって計算できる
> AIC(mr)
[1] 638.9382
> BIC(mr)
[1] 648.189
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
28 / 56
複数のモデルの推定
方法 1 ：今までの作業をくりかえす
ひとつひとつモデルを手作業で推定していく
このような方法を対話式という
方法 2 ：バッチファイルを用いる
必要なコマンドをファイルに書きこむ
結果をまとめて最後に出力するコマンドも
ファイルに書いておく
ファイルを実行する
このようなファイルをバッチファイルという
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
29 / 56
バッチファイルによる方法
ファイル → 新しいスクリプト
エディターに次のように書きこむ
cvdata <- read.csv("D:/Rdata/cv.csv")
model1 <- lm(conv ˜ pop + car + nm_f_30, data=cvdata)
aic1 <- AIC(model1)
model2 <- lm(conv ˜ pop + gpp + divorce, data=cvdata)
aic2 <- AIC(model2)
model3 <- lm(conv ˜ pop + sm + w_pop, data=cvdata)
aic3 <- AIC(model3)
print(c(aic1,aic2,aic3))
適当なファイル名で適当な場所に保存
編集 → 全て実行
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
30 / 56
演習 2：コンビニ件数のモデル
演習 2
すべての変数を用いてコンビニ件数のモデルとし
て, AIC, あるいは BIC の意味で最適なものを探索
し , 結果選択されたモデルについて考察せよ.
.
H. Hara (Niigata U.)
R による計量経済分析
.Nov 11-, 2014
31 / 56
Outline
1
都道府県別コンビニ数のデータ
記述統計
回帰分析
2
物価と為替レートの関係
3
賃金関数の推定
4
たばこの価格と消費量のデータ
5
レポート課題
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
32 / 56
データのダウンロード
講義のサイトからデータ”cpi.csv” をダウンロード
して Rdata に保存
http://www.econ.niigata-u.ac.jp/˜hara/grad-ecm/
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
33 / 56
データと分析の目的
1973 年∼2013 年までの日本の
.
.
.
1
2
.
3
CPI : 消費者物価指数
EXP : 円 /ドルレート
WAGE : 賃金指数
データ元
.
.
1
消費者物価指数：総務省統計局
http://www.stat.go.jp/data/cpi/
2
円 /ドルレート：日銀時系列データ検索サイト
http://www.stat-search.boj.or.jp/
.
.
3
賃金指数 : 毎月勤労統計調査 (厚労省)
http://www.mhlw.go.jp/toukei/list/30-1.html
長期時系列表
就業形態別現金給与総額就業形態計 (30 人以上)
.
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
34 / 56
データの読み込み
データを cpidata に読み込む
> cpidata <- read.csv("D:/Rdata/cpi.csv")
各変数をオブジェクトに入力
> cpi <- cpidata$CPI
> exr <- cpidata$EXR
> wage <- cpidata$WAGE
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
35 / 56
物価と為替レートの関係
通常の経済理論によれば・
・
・
円高 (ドル安)
⇒ 輸入品物価の下落
⇒ 国産品の物価も下落
円高 (ドル安) ⇔ 円 /ドルレートは小
円 /ドルレートと消費者物価指数は比例的であること
が予想できる
cpii : 消費者物価指数 (CPI), exri : 円 /ドルレート
cpii = α + βexpi + ui
という単回帰モデルを推定すると , β > 0 となるはず .
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
36 / 56
単回帰分析
OLS で推定
> cpi.sr <- lm(cpi ˜ exr)
> result.sr <- summary(cpi.sr)
散布図と推定直線を描く
> plot(cpi ˜ exr)
> abline(cpi.sr,col=’’red’’)
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
37 / 56
推定結果
βˆOLS < 0 で有意
理論が誤りなのか？
cpi
40
50
60
70
80
90
100
Estimate Std. Error t value Pr(>|t|)
(Intercept) 123.89027
2.87927
43.03 < 2e-16 ***
exr
-0.21354
0.01672 -12.77 1.65e-15 ***
100
150
200
250
300
exr
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
38 / 56
物価と為替レートの関係
物価は為替レートだけでなく他の要因にも依存
賃金指数 (単位時間当の賃金を指数化したもの)
賃金率高 → 物価上昇
WAGEi : i 年の賃金指数
cpii = α + β2 exri + β3 wagei + ui
を推定してみよ
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
39 / 56
重回帰分析
OLS で推定
> cpi.mr <- lm(cpi ˜ exr + wage)
> result.mr <- summary(cpi.mr)
> result.cpi$coefficients
推定結果
Estimate Std. Error
t value
Pr(>|t|)
(Intercept) 18.12786174 5.99872038 3.0219548 4.478012e-03
exr
0.00919046 0.01364285 0.6736465 5.046132e-01
wage
0.78348265 0.04387526 17.8570502 4.471507e-20
βˆ2 > 0
単回帰モデルは過少定式化だったことが疑われる
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
40 / 56
演習 3：モデルの評価
演習 3
推定した 2 つのモデルの AIC, BIC を計算して, そ
の結果について考察せよ.
.
.
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
41 / 56
Outline
1
都道府県別コンビニ数のデータ
記述統計
回帰分析
2
物価と為替レートの関係
3
賃金関数の推定
4
たばこの価格と消費量のデータ
5
レポート課題
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
42 / 56
データのロード
R には多数のサンプルデータが用意されている.
そのリストは
> data()
で見ることができる.
AER の CPS1988 をロード
> data(CPS1988)
> CPS1988
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
43 / 56
CPS1988 データ
Current population survey 1988 in US
18 歳∼70 歳で $5 以上の収入のある 28,155 人の
男性の
wage : 賃金
education : 教育年数
experience : 就業年数
ethnicity : 人種 (白人 or 黒人)
smsa : 居住地が大都市統計地域か否か
region : 居住地域 (北西部, 中西部, 南部 or 西部)
parttime : 常勤 or 非常勤
下の 4 つは質的変数
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
44 / 56
賃金関数の推定
以下のモデルを推定してみよう.
log(wage) = α + β2 experience
+ β3 experience2 + β4 education
+ β5 ethnicity + ui
> cps_lm <- lm(log(wage) ˜ experience + I(experienceˆ2)
+ education + ethnicity, data=CPS1988)
>
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
45 / 56
推定結果
> result.cps_lm <- summary(cps_lm)
> result.cps_lm$coefficients
Estimate
Std. Error
t value
Pr(>|t|)
(Intercept)
4.321394996 1.917421e-02 225.37534 0.000000e+00
experience
0.077473231 8.800466e-04 88.03310 0.000000e+00
I(experienceˆ2) -0.001316066 1.898751e-05 -69.31223 0.000000e+00
education
0.085672819 1.272186e-03 67.34298 0.000000e+00
ethnicityafam
-0.243364296 1.291812e-02 -18.83898 1.104188e-78
>
ethnicity は
{
ethnicityafam =
1, 黒人
0, 白人
というダミー変数に変換される.
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
46 / 56
賃金関数の推定
以下のモデルを推定してみよう.
log(wage) = α + β2 experience
+ β4 education + ui
> cps_lm2 <- lm(log(wage) ˜ experience + education,
data=CPS1988)
>
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
47 / 56
anova を用いた線形制約の検定
H0 : β3 = β5 = 0 の検定
> anova(cps_lm2,cps_lm)
Analysis of Variance Table
Model 1: log(wage) ˜ experience + education
Model 2: log(wage) ˜ experience + I(experienceˆ2)
+ education + ethnicity
Res.Df
RSS Df Sum of Sq
F
Pr(>F)
1 28152 11357.5
2 28150 9598.6 2
1758.9 2579.2 < 2.2e-16 ***
--Signif. codes: 0 ‘ *** ’ 0.001 ‘ ** ’ 0.01 ‘ * ’ 0.05 ‘ . ’ 0.1 ‘
1
>
’
RRSS = 11357.5, U RSS = 9598.6
RRSS − U RSS = 1758.9
F = 2579.2
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
48 / 56
演習 4：賃金関数の推定
演習 4
他の変数も用いて, 賃金関数として適切なモデル
を探索し , 結果について考察せよ.
.
region を用いると, 3 つ以上のグループを持つ
ダミー変数の変換のされ方がわかる.
.
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
49 / 56
Outline
1
都道府県別コンビニ数のデータ
記述統計
回帰分析
2
物価と為替レートの関係
3
賃金関数の推定
4
たばこの価格と消費量のデータ
5
レポート課題
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
50 / 56
たばこ消費の価格弾性値の推定
講義のサイトからデータ”US-CigarCons.csv” を
DL して, D:/Rdata/に保存
http://www.econ.niigata-u.ac.jp/˜hara/grad-ecm/
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
51 / 56
たばこ消費の価格弾性値の推定
1985 年から 1995 年までのアメリカの 48 州 (アラ
スカ, ハワイ以外) のたばこの価格と消費量に関す
るデータ
アメリカは州ごとにたばこ税も消費税も異なる
たばこ消費量の価格弾性値が推定できる
データの出典は Stock and Waton の計量経済の
教科書のサポートサイト
http://wps.aw.com/aw_stock_ie_3/178/45691/11696965.
cw/index.html
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
52 / 56
データに含まれる変数
state : 州
year : 年 (1985∼1995)
cpi : 消費者物価指数
pop : 州の人口
packpc : 一人当の消費量 (パック数)
income : 州の総個人所得
tax : excise tax(たばこ税, $ 0.01)
avgprs : たばこ 1 パックの平均価格
(含消費税, $ 0.01)
taxs : 消費税 (含たばこ税, $ 0.01)
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
53 / 56
演習 5：たばこ消費の価格弾性値の推定
演習 5
1
弾性値モデル
log packpci = α + β log avgprs + ui
2
を推定して, 結果について考察せよ.
他に説明変数を加えて, AIC, BIC の意味でより好
ましいモデルを探索し , 結果選ばれたモデルに基
づいて, たばこ消費と価格の関係について考察
.
せよ.
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
54 / 56
Outline
1
都道府県別コンビニ数のデータ
記述統計
回帰分析
2
物価と為替レートの関係
3
賃金関数の推定
4
たばこの価格と消費量のデータ
5
レポート課題
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
55 / 56
課題
課題
Baltagi の教科書 p.90 の問 15 を R を用いて解答
せよ
データは講義のサイトの usgas.csv を用いよ.
Baltagi の教科書 p.90 の問 16 を R を用いて解答
せよ.
データは講義のサイトの natural.csv を用いよ.
.
H. Hara (Niigata U.)
R による計量経済分析
Nov 11-, 2014
56 / 56

Download Report