Rによる単回帰分析

Rによる単回帰分析
自己回帰モデルへの橋渡し
高崎経済大学宮田庸一
Rとは
• Rとは、統計処理、グラフ作成のために開発され
た言語である。
• Rはフリーソフト
• ダウンロードおよびインストールの仕方
→高崎経済大学宮田研究室のRというページ
• 高崎経済大学のRは現在のところ,
英語版しか入っていない・・・(-_-;)
データのインポート
• GDP.csvをRに取り込む.
• Rにデータを取り込むために
data1<-read.csv("C:\\Import\\GDP.csv",header=T,row.names=1)
と入力しEnterを押す
パスを指定する
1980～2009までの実質民間
最終消費支出（Private Final
Consumption Expenditure）
と実質GDP:
http://www.esri.cao.go.jp/jp/s
na/qe1022/gdemenu_ja.html
3
GDP（10億円)
PFCE （10億円)
1980
287366.4
167753.7
1981
298687.1
171754.8
1982
308057
179712.9
1983
318921.7
185174.5
1984
334110.7
190722.2
1985
355096.2
199114.5
1986
361807.1
206209.7
データのグラフ表示
散布図：plot(data1,main="消費と支出")
ヒストグラム:hist(data1$GDP,main="GDP")
GDP
260000
220000
民間最終消費支出(PFCE)
6
4
180000
2
0
Frequency
8
300000
10
消費と支出
250000
300000
350000
400000
450000
data1$GDP
500000
550000
600000
300000
350000
400000
450000
GDP
500000
550000
データを指定
単回帰分析
回帰分析の実行
回帰分析の出力
result1<-lm(PFCE~GDP,data=data1)
summary(result1)
-----Rでの出力例-----------------------------------------------> summary(result1)
Call:
lm(formula = PFCE ~ GDP, data = data1)
Residuals:
Min
1Q Median
3Q Max
-6695.1 -2768.7 196.3 2445.7 8813.3
y  8.126103  0.5465x
標準誤差 t値
P値
Coefficients: 推定値
Estimate
Std. Error
t value Pr(>|t|)
(Intercept) 8.126e+03 3.941e+03 2.062
0.0486 *
GDP
5.465e-01 8.585e-03 63.654 <2e-16 ***
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
決定係数
Residual standard error: 3834 on 28 degrees of freedom
Multiple R-Squared: 0.9931, Adjusted R-squared: 0.9929
F-statistic: 4052 on 1 and 28 DF, p-value: < 2.2e-16
自由度調整済み決定係数
回帰直線の上書き
abline(result1)
他にもモデル診断もあるが省略する
260000
220000
180000
民間最終消費支出(PFCE)
300000
消費と支出
300000
350000
400000
450000
GDP
500000
550000
予測
200000
result.pre
直線に見えるが、実は曲線
250000
300000
• GDPが400000の時, 民間最終消費支出(PFCE)を予測する.
> new<-data.frame(GDP=400000)
> predict(result1,new,interval="predict",level=0.95)
fit
lwr
upr
回帰直線と区間予測
[1,] 226718.4 218684.2 234752.6
300000
350000
400000
450000
new$GDP
500000
550000
時系列データ
• 為替、株価、マクロ統計、心電図、筋電図etc
時刻tによりランダムに変動するデータを時系
列データと言う.
, X 1, X 0 , X1, X 2 ,....,X T ,....
• 時系列データを通常の統計手法で処理するた
めには，定常性を仮定しなくてはいけない．
• 定常性とは,
E( X t )  , Cov( X t , X t s )   ( s)
時系列分析の目的
• どこまで、過去のデータに影響を受けるのか
を推定
• 株式データ、為替データの場合、リスクがど
のくらいあるのかの推定
• 過去のデータから将来の値を予測する
時系列データのプロット
94
2010/1/4から2010/10/25までの円ドル為替相場の週足データ(rate_data.csv)
data2<-read.csv("C:\\Import\\rate_data.csv",header=T,row.names=1)
plot(data2$doll.yen,type="l",xlab="時刻",ylab="為替レート")
88
86
84
82
80
為替レート
90
92
これは定常ではないと考え
られる.
0
10
20
時刻
30
40
データの差分を取る
data3<-diff(log(data2$doll.yen))
plot(data3,type="l“,xlab="時刻",ylab="差分したデータ")
定常でないデータ
0.02
, X 1, X 0 , X1, X 2 ,....,X T ,....
0.00
-0.01
-0.02
,Y1,Y0 ,Y1,Y2 ,....,YT ,....
-0.03
定常なデータとなる
差分したデータ
Xt
Yt  log
 log X t  log X t 1
X t 1
0.01
差分を取る
0
10
20
時刻
30
40
自己相関係数(acf)
• 時系列データ
y1 , y2 ,..., yT
1 T
1 T
y   yt , ˆ (0)   ( yt  y ) 2
T t 1
T t 1
T
1
ˆ ( s)   ( yt  y )( yt  s  y )
T t  s 1
t
ドル/円
ドル/円
t-1
1月11日
90.7799
92.61
1月4日
1月18日
89.8899
90.7799
1月11日
1月25日
90.3
89.8899
1月18日
2月1日
89.33
90.3
1月25日
2月8日
90.0199
89.33
2月1日
2月15日
91.5899
90.0199
2月8日
2月22日
88.8499
91.5899
2月15日
3月1日
90.29
88.8499
2月22日
3月8日
90.4599
90.29
3月1日
• 自己相関係数(Auto-correlation)
T
(y
ˆ ( s ) t  s 1
ˆ ( s ) 

ˆ (0)
t
 y ) ( yt  s  y )
T
(y
t 1
t
 y)
2
Rでの実行
ˆ (0), ˆ (1),...,ˆ (20) をプロット
acf(data3, main="ACF")
0.4
0.2
0.0
-0.2
ACF
0.6
0.8
1.0
ACF
0
5
10
Lag
15
1次の自己回帰モデル
Y1 , Y2 ,...,YT
Yt  Yt 1  ut
(t  1,...,T )
ここをXtと考えると, 単回帰モデルと同じように推定が出来
そう.
→Xtは定数ではなく，確率変数であるが，実はこの推定法
はある条件の下で正しいことが知られている．
2
• パラメーター
  ( , )
推定法
• Yule Walker法（山本 p57）
T
ˆ 
 (Y  Y )(Y
t 2
t
t 1
Y )
T
2
(
Y

Y
)
 t
t 1
T
1 T
1
2
ˆ   (Yt  Y )  ˆ  (Yt  Y )(Yt 1  Y )
T t 1
T t 2
2
1 T
Y  Yt
T t 1
Rでの実行(ar)
ar01<-ar(data3,aic=FALSE,order.max=1)
> ar01
Call:ar(x = data3, aic = FALSE, order.max = 1)
Coefficients:
1
-0.2281
Yt  0.2281Yt 1  ut
ut～N (0,0.0001944)
Order selected 1 sigma^2 estimated as 0.0001944
金融工学においては, 0.0001944をリスクと考える.
1次の自己回帰モデル
• 切片のある時系列モデルも考えられる
Yt    Yt 1  ut
(t  1,...,T )
ここをXtと考えると, 単回帰モデルと同じように推定が出来
そう.
2
• パラメーター
  ( ,  , )
Rでの実行
data2<-read.csv("C:\\Import\\rate_data.csv",header=T,row.names=1)
data3<-diff(log(data2$doll.yen))
ar1<-arima(data3, order=c(1,0,0))
ARIMA Modelという，より一般的なモデルがある.
> ar1
Call:
arima(x = data3, order = c(1, 0, 0))
Yt  0.0033 0.2328Yt 1  ut
ut～N (0,0.0001847)
Coefficients:
ar1 intercept
-0.2328 -0.0033
s.e. 0.1514 0.0017
sigma^2 estimated as 0.0001847: log likelihood = 120.91, aic = -235.83
参考文献
[1]平成20年度国民経済計算（平成12年基準・
93SNA）http://www.esri.cao.go.jp/jp/sna/h20kaku/22annual-report-j.html
[2] 山本拓, (2009), 経済の時系列分析，創文
社
[3] 船尾暢男, (2009), The R tips, オーム社

Download Report