見せかけの回帰と共和分回帰

options(width=100,warn=-1,digits=5,messages=1)


#計量経済学のテキストでは必ず取り上げられているテーマであるが
#ランダムウォーク過程に従う仮想データの生成法は分析の目的により
#よって様々な方法があるようだ。任意の初期値がその時系列にどのよう
#な影響を見たい場合もあり、あるいは、初期値の影響度が薄まった時期
#のデータを使ってシミュレーションを行うこともある。例えば、実際
#の計算に必要なデータが１０００個としたときに５０個大きめにして
#１０５０個の疑似データを生成させてから、最初の５０個を除外した
#残りの１０００個の生成データを使ってシミュレーションをする仕方で
#ある。実証のための計量時系列分析　エンダース著　新谷、藪　訳
#有斐閣　２０１９ではブートストラップ法の解説で、５０という数値
#にこだわらずに初期値の影響をみて決めればよいと述べられている。
#こでは幾つかの簡単なランダムウォーク過程の生成法をみてみる。
#ライブラリーとしては以下のものを使うことにする。

suppressPackageStartupMessages(library(lmtest))

#ダービンワトソン統計量を計算するのに使う
suppressPackageStartupMessages(library(tseries))

#単位根テストで使う

##単純なランダムウォーク系列の生成

set.seed(10)
t<-c(1:300) 

# 300個の正規乱数（平均値０、標準偏差　１）を作る

x<-ts(rnorm(n=length(t)))
## cumsum関数を使ってランダムウォーク過程を生成する
xw<-cumsum(x)
ts.plot(xw,type="l")






# forを使った繰り返し文で同様のデータを作成することもできる
##############
set.seed(10)
t<-c(1:300)
y <- w <- ts(rnorm(length(t)))
for (t in 2:length(t) ) y[t] <- y[t-1] + w[t]
ts.plot(y, type="l")





head(xw)

## [1]  0.018746 -0.165506 -1.536837 -2.136005 -1.841460 -1.451665head(y)
## [1]  0.018746 -0.165506 -1.536837 -2.136005 -1.841460 -1.451665
#いずれも同一のランダムウォーク系列となっている。

#たとえば初期値を０から始めるランダムウォーク系列を作る場合には
##cumsum関数を使う

t<-c(0:300)

# 300個の正規乱数（平均値０、標準偏差　１）
set.seed(10)
x1<-ts(rnorm(n=length(t)-1))
x1w<-c(0,cumsum(x1))
plot(x1w,type="l")





#forを使った繰り返し文

set.seed(10)
t<-c(0:300)
x2<-w2<-rep(0,length(t))
x2 <- w2 <-ts(c(0, rnorm(n=length(t)-1)))
for (t in 2:length(t) ) x2[t] <- x2[t-1] + w2[t]
ts.plot(x2, type="l")




head(x1w)
## [1]  0.000000  0.018746 -0.165506 -1.536837 -2.136005 -1.841460head(x2)
## [1]  0.000000  0.018746 -0.165506 -1.536837 -2.136005 -1.841460
# いずれも初期値０とする同一の系列となっている


#見せかけの回帰について
#最初にx,yのランダムウォーク過程のデータを作ると
set.seed(10)
x<-ts(cumsum(rnorm(n=300)))
set.seed(30)
y<-ts(cumsum(rnorm(300)))

#次に、ランダムウォーク過程に従うｘとｙについて単純回帰分析をする
lmmod1 <- lm(y ~ x)
summary(lmmod1)

##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -20.635  -6.309  -0.561   6.348  21.596
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   8.9761     1.3985    6.42  5.4e-10 ***
## x             1.8643     0.0772   24.16  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.54 on 298 degrees of freedom
## Multiple R-squared:  0.662,  Adjusted R-squared:  0.661
## F-statistic:  584 on 1 and 298 DF,  p-value: <2e-16

dwtest(lmmod1)

##
##  Durbin-Watson test
##
## data:  lmmod1
## DW = 0.0602, p-value <2e-16
## alternative hypothesis: true autocorrelation is greater than 0

#ｘとｙは全く相互に独立している無関係のデータであるが、決定係数
#は0.6621と高い値となっている。ｘ、ｙ共に非定常過程であるランダム
#ウォーク過程に従っているが、両者の確率的トレンドに強い相関が
#あったため見せかけの回帰現象が生じている。
#Damodar GujaratiのBasic Econometrics からの受け売りとなってしまう
#がGranger-Newboldによれば決定係数がダービンワトソン統計量よりも大
#きい時には経験則として見せかけの回帰を疑うことができるとされてい
#る。
#この事例でも決定係数は高い割にはDWは０．０６と低い。
#見せかけの回帰が疑われるときには、各変数の一次差分を計算して
#それについて回帰分析をすることで識別できる可能性がある。

resid<-ts(lmmod1$residuals)

plot(resid,type="l")






acf(resid)




#上記のように残差系列residをグラフで見るとランダムウォーク
#過程のように見える。また残差のコレログラムを見ると系列相関
#係数が#なだらかに減少しており、誤差の影響が長期にわたって
#記憶されることを示している。
#この残差系列が定常かどうかをadf.testで単位根検定をしてみる。なお、
#単位根検定の考え方については別途にRを使ってディッキー・フラー
#DF分布のシミュレーションを行ってみた。（単位根検定、
#　Rでディッキー・フラーDF分布をシミュレート
#　をご参照。

adf.test(resid

)##
##  Augmented Dickey-Fuller Test
##
## data:  resid
## Dickey-Fuller = -2.74, Lag order = 6, p-value = 0.27
## alternative hypothesis: stationary

# ライブラリー tseries のマニュアルによればadf.testでは帰無仮説
#は単位根あり、つまり非定常過程であると設定されている。
#仮に有意水準を５％と置いたとすれば、p値が0.266と有意水準の
#比べて大きいので、帰無仮説である非定常変数を棄却できない。
#従って、残差系列residは非定常過程に従う#のだろうと暫定的に
#判断することになる。

#ｘとｙの一時差分について回帰分析をすると

lmmod2 <- lm(diff(y) ~ diff(x))
summary(lmmod2)

##
## Call:
## lm(formula = diff(y) ~ diff(x))
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -2.8041 -0.7417 -0.0825  0.7231  2.8116
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|) 
## (Intercept)  -0.1393     0.0617   -2.26    0.025 *
## diff(x)      -0.0138     0.0640   -0.22    0.829 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.06 on 297 degrees of freedom
## Multiple R-squared:  0.000157,   Adjusted R-squared:  -0.00321
## F-statistic: 0.0466 on 1 and 297 DF,  p-value: 0.829

dwtest(lmmod2)

##
##  Durbin-Watson test
##
## data:  lmmod2
## DW = 2.1, p-value = 0.8
## alternative hypothesis: true autocorrelation is greater than 0

#今度は決定係数はゼロに近い値となっている。ダービンワトソン統計量
#も２に近く系列相関の低いことを示している。つまりｘとｙは相互の
#関係性が低いと判断される。

#xとyの一次差分について単位根テストをしてみると
adf.test(diff(x))

##
##  Augmented Dickey-Fuller Test
##
## data:  diff(x)
## Dickey-Fuller = -6.5, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationaryadf.test(diff(y))
##
##  Augmented Dickey-Fuller Test
##
## data:  diff(y)
## Dickey-Fuller = -6.72, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
#　p値が0.01より小さく有意水準を大きく下回っており帰無仮説の
#非定常変数を棄却できる。つまりｘの一次差分、yの一次差分は
#定常過程に従うと判断される。


#共和分回帰について

#次にランダムウォ－ク過程に従うデータ、ミュー(mu)を作り、
#そのmuに互いに独立なホワイトノイズを加えたx1,y1の２つの
#非定常系列を生成する。

T<-300
mu<-ts(cumsum(rnorm(T)))
plot(mu)





x1<-ts(mu+rnorm(300))
y1<-ts(mu+rnorm(300))

#x1とy1はランダムウォーク過程に従うので非定常となるが
#両者ともに非定常変数ｍｕが含まれている。
#下記のグラフからy1 x1ともに似たようなランダムウォーク
#過程を示している
matplot(cbind(x1,y1),type="l",lwd=2)





#y1のx1への回帰直線を求めると
lmmod3 <- lm(y1 ~ x1)
summary(lmmod3)

##
## Call:
## lm(formula = y1 ~ x1)
##
## Residuals:
##    Min     1Q Median     3Q    Max
## -4.762 -1.047  0.017  0.986  3.913
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   0.0857     0.0869    0.99     0.33   
## x1            0.9503     0.0152   62.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.47 on 298 degrees of freedom
## Multiple R-squared:  0.929,  Adjusted R-squared:  0.929
## F-statistic: 3.93e+03 on 1 and 298 DF,  p-value: <2e-16

#決定係数も高く係数の推定値も有意となっているのが
#これは見せかけの回帰になるのだろうか。
#x1 y1　それぞれについてADF　テストをすると下記のように
#x1 y1 は非定常であることを棄却できない、つまり暫定的に
#非定常だろうと判断される。

adf.test(x1)

##
##  Augmented Dickey-Fuller Test
##
## data:  x1
## Dickey-Fuller = -1.18, Lag order = 6, p-value = 0.91
## alternative hypothesis: stationaryadf.test(y1)##
##  Augmented Dickey-Fuller Test
##
## data:  y1
## Dickey-Fuller = -0.856, Lag order = 6, p-value = 0.96
## alternative hypothesis: stationary

#x1とy1が共和分の関係なるかどうか、Phillips-Ouliaris
#テスト(po.test）をする。

po.test(cbind(x1,y1))

##
##  Phillips-Ouliaris Cointegration Test
##
## data:  cbind(x1, y1)
## Phillips-Ouliaris demeaned = -307, Truncation lag parameter = 2, #p-value = 0.01
# ライブラリー tseries の　po.testでは帰無仮説は共和分で
#はないと設定されている。このテスト結果ではｐ値が1％以下と
#なっており、共和分でないとする帰無仮説は棄却される。
#つまり、x1 y1 は共和分の関係にある判断できる。
#このため、上記の回帰分析はx1を説明変数、y1を従属変数と
#する共和分回帰を行ったことになる。
#共和分（cointegration）については、どの計量経済学のテキスト
#でも詳しく論じられているが、大ざっぱに言えばn個の系列、例えば
#ｘ、ｙ、の二つの#系列の線形結合が定常過程となることをいう。
#先ほどの見せかけの回帰ｘ、ｙの例では
#ｘとｙの線形結合y-8.96-1.86x=residは非定常であったが
#x1とy1の線形結合y1-0.085-0.95x1=resid3は定常過程に従っている。
#この共和分回帰の残差系列resid3について、念のためにADFテストを
#してみる。下記のようｐ値が１％で非定常とする帰無仮説は
#棄却される。つまり、残差系列は定常過程に従っており、x1とy1の
#線形回帰分析は意味があることになる。
#残差系列resid3のコレログラムを見ても系列相関係数は急激に低下し
#誤差の影響は長期に及ばないことがわかる。

resid3<-ts(lmmod3$residuals)
adf.test(resid3)

##
##  Augmented Dickey-Fuller Test
##
## data:  resid3
## Dickey-Fuller = -5.6, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationaryacf(resid3)





#このように共和分回帰の場合は見せかけの回帰と異なり、x1とy1
#の間に何らかの長期的な均衡関係が存在することを示唆している。
#主要国のGDPなどのマクロ経済指標、主要な穀物等の商品価格、
#各種の金利、為替相場等の時系列データは共和分関係が見られる
#ことがあるようだ。
#もし、うまい具合に株価指数や為替レートについて共和分回帰で
#長期的均衡を見つけることが出来れば,目先の短期的な乖離は均衡線
#に戻ることを期待する投資戦略も考え得るが,成功するかどうかは
#わからない。均衡値に戻るのに余りにも時間がかかれば実用的でないだろう。