単位根検定の背景（ディッキー・フラー分布のシミュレーション with R）

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

#ディッキー(Dickey)とフラー(Fuller)はランダムウォーク過程では最小
#２乗推定量にバイアスが生じることを示した。
#単位根検定（unit root test）(ADF検定）ディッキーフラー(Dickey-Fuller)
#検定の考え方と背景について簡単なシミュレーションとグラフで検討する。
#まずはじめに、ランダムウォーク過程の分散について見てみる。
#ランダムウォーク過程の分散
#ランダムウォーク過程に従う時系列データの特徴をグラフで
#見てみる。u(t)は期待値０、分散σ²　のホワイトノイズとする。
#ホワイトノイズは系列相関がゼロを前提にするが、i.i.d(independently
# and identically distributed)のようにu(t)が独立であることまで
# 要求しないという意味で、i.i.dより弱い条件となっている。（相関が
#ゼロであったとして必ずしも独立とは限らないという意味で）
#ランダムウォーク過程を　ｙ（ｔ）＝ｙ（ｔ－１）＋u(t)　で表すと　
#  E(u(t))=0   E(u(t)²)=σ²  E(u(t)・u(t-i))=0 i≠0 の前提から
# Var(y(t))=t×σ²　E(y(t))=0 となることが、どの時系列分析のテキスト
#にも書かれている。ランダムウォーク過程に従う時系列はその分散が一定
#でなく時間の経過と共に増加していく特徴がある。
#この現象を正規分布（平均値ゼロ、標準偏差　１）に従うホワイトノイズ
#によるシミュレーションで表すことを試みる。
#手順としては１００期間（例えば時間であ１００時間）の時系列を
#２０００本作成してその時間の経過後との平均値と分散の変化を視覚化して
#みる。
#

rw <- list()  # 乱数のリストを準備
set.seed(1)

Nt<-100  # サンプル数（期間数）

M<-2000 #繰返し回数
# 平均値０　標準偏差 １　の正規分布乱数を
#100期間の時系列を2000本作成する

for(i in 1:M){rw[[i]] <- rnorm(Nt)}

# リストを行列に変換する
rwmat<-matrix(unlist(rw),nrow=Nt,byrow=F)
#行列に初期値０の行を追加する
rwalk<-rbind(rep(0,M),rwmat)


#２０００本のランダムウォーク系列を作成する
rwalk1<-apply(rwalk,2,cumsum)
#２０００本のランダムウォーク系列のグラフ
matplot(rwalk1,type="l",xlab="TIME")







# 経過期間ごとの平均値を計算する
rwalkmu<-apply(rwalk1,1,mean)
ts.plot(rwalkmu,type="l",ylim=c(-1,1))






#大体ゼロ近辺で推移している。シミュレーションの反復計算
#回数を1万とかに増やせば直線に近くなるがPCの計算時間がかかる。
# 経過期間ごとの分散を計算する
rwalkvar<-apply(rwalk1,1,var)
ts.plot(rwalkvar,type="l")




#大体、時間の経過と共に分散が増加する傾向で推移している。シミュレーション
#の反復計算回数を1万とかに増やせば直線に近くなるがPCの計算時間がかかる。
#
#最小2乗推定量のバイアスについて
#このように分散が時間の経過と共に増加する時系列データに
#最小二乗法を適用してパラメータを推定した場合にはバイアスが生じる
#ことが知られている。回帰分析などを行った場合には時間の経過と共に
#分散が変化するのでｔ分布の形も変わってしまう。
#ディッキー(Dickey)とフラー(Fuller)はランダムウォーク過程では最小２乗
#推定量にバイアスが生じることを示した。
#真のモデルがｙ（ｔ）＝ρ・ｙ（ｔ－１）とした場合に
#両辺からｙ（ｔ－１）を差し引いて
#ｙ（ｔ）－ｙ（ｔ－１）＝△ｙ（ｔ）,β＝ρ－１　とすれば
# △ｙ（ｔ）＝β・ｙ（ｔ）　で表現できる。この式について
#最小二乗法を適用してβを推定するとρが1であれば、つまり単位根が存在
#すればβはゼロとなることが予想される。
#
#そこで3種類のランダムウォーク
#（定数項（ドリフト　α0）が無い場合
#△ｙ（ｔ）＝β・ｙ（ｔ）+ u(t)
#ドリフト　α0　がある場合
#△ｙ（ｔ）＝α0 + β・ｙ（ｔ）+ u(t)
#、トレンド項ｔとドリフトα0の両者がある場合
#△ｙ（ｔ）＝α0 + β・ｙ（ｔ）+ α１・ｔ＋u(t)
#の三種類のケースについて最小２乗法を使って回帰係数とｔ値を推定してみる。

########################
#ディッキー、フラー（DF）分布のシミュレーション
#による近似計算
nt<-100  #サンプル数
nfirst = 200　　#多めにとって、その後の回帰分析では使用しないサンプル数
nsmp = nt + nfirst + 1  #シミュレーションに必要なサンプル数
nrep = 1000   # シミュレーションの反復計算回数
#３種類のモデルの回帰係数をプールするベクトルを
# bvalnoconst,bvaldrift,bvaltrend で表す
#
#同様にして回帰係数の推定値に関するｔ値をプールするベクトルを
#tvalnoconst、-tvaldrift、tvaltrend　で表す。

bvalnoconst<-bvaldrift<-bvaltrend<-rep(0,nt) 
tvalnoconst<-tvaldrift<-tvaltrend<-rep(0,nt)　
#下記でｙはランダムウォーク系列y(t)を示し
#yyはy(t-1)　つまり1期遅れの系列を示し
#dyは△ｙ（ｔ）を示している。
#あとは回帰分析の結果（summary)から必要な数値を抽出する
y<-YY<-dy<-rep(0,nrep)
set.seed(1)
for (i in 1:nrep) {
  y <- cumsum(rnorm(nsmp))
  yy<-y[-c(1:nfirst)]
  dy <- diff(yy)
 
  lmnoconst <- summary(lm(dy ~ 0 + yy[-length(yy)]))
  tvalnoconst[i]<- lmnoconst$coefficients[1,3]
  bvalnoconst[i]<- lmnoconst$coefficients[1,1]
  lmdrift <- summary(lm(dy ~ yy[-length(yy)]))
  tvaldrift[i] <- lmdrift$coefficients[2,3]
  bvaldrift[i] <- lmdrift$coefficients[2,1]
  lmtrend <- summary(lm(dy ~yy[-length(yy)] + c(1:nt)))
  tvaltrend[i] <- lmtrend$coefficients[2,3]
  bvaltrend[i] <- lmtrend$coefficients[2,1]
}

#ドリフトなし のモデルのβの分布

plot(density(bvalnoconst),col="blue",lwd=5,main="ドリフトなし　βの分布")




#真のモデルがランダムウォーク過程に従っていれば、ρ＝１つまり単位根が
# 存在するのでβ＝ρ－１＝０　となる筈であるがβの推定値はゼロより小さい
#マイナス側に偏っていることが分かる。ドリフト付き、トレンド付きになると
#マイナス側への偏りは、さらに大きくなる。仮にβ＝－０．１となれば
#ρ＝０．９であり、単位根が真に存在するときでも最小２乗推定値は０．９
#と計算する可能性を意味している。

plot(density(bvaldrift),col="blue",lwd=5,main="ドリフト付き　βの分布")



plot(density(bvaltrend),col="blue",lwd=5,main="ドリフト、トレンド付き　βの分布")


#３種類もモデルについてのｔ値の分布をグラフにすると左右対照のｔ分布
#とは異なっておりマイナスである側に偏っていることがわかる。
#単位根の検定をするときには通常のｔ分布表を使うことは適当でなく
#マイナス側への偏りを考慮する必要があることがわかる。ざっくり言えば
#より大きいマイナスのｔ値が出てくれば帰無仮説の単位根あり（非定常）
#を棄却できることがわかる。


plot(density(tvalnoconst),col="green",lwd=5,main="ドリフトなし　ｔ値の分布")




plot(density(tvaldrift),col="green",lwd=5,main="ドリフト付き　ｔ値の分布")





plot(density(tvaltrend),col="green",lwd=5,main="ドリフト、トレンド付き　ｔ値の分布")




print("ドリフトなし　　サンプル数　１００")
## [1] "ドリフトなし　　サンプル数　１００"
quantile(tvalnoconst,c(0.01,0.05,0.1))
##      1%      5%     10%
## -2.5517 -1.9384 -1.6248
print("ドリフト付き　　サンプル数　１００")
## [1] "ドリフト付き　　サンプル数　１００"
quantile(tvaldrift,c(0.01,0.05,0.1))
##      1%      5%     10%
## -3.6705 -2.9017 -2.6580
print("ドリフト、トレンド　付き　　サンプル数　１００")
## [1] "ドリフト、トレンド　付き　　サンプル数　１００"
quantile(tvaltrend,c(0.01,0.05,0.1))
##      1%      5%     10%
## -4.0456 -3.4542 -3.1615
#βの推定値の分布はゼロより少し偏小さい数値に偏って分布している。
# 非定常過程であればρ＝１なのでｂ＝０　を中心にした分布が
#期待されるが、そうなっていない。βに係るｔ値の分布を見てみると
#左右対象のｔ分布ではなくゼロより小さい方に偏っている。
#3種類のモデルについてのｔ値の分位表を見ると単純なシミュレーション
#による結果だが、ディッキーフラーのｔ検定でのタウ統計量τ=100の場合
#のDF検定表の数値に近いものとなっている。
#時間がかかるが反復計算回数を５万ぐらいにするともっとDF検定表に
#近づく。しかし、実際のところは、Rであれば様々なパッケージで
#単位根検定(ADF検定）の関数が用意されているので、特に自らが検定表を作
#る必要はないが単位根検定(ADF検定）の関数が大ざっぱにどのようなことを
#しているのかを知っておくのも損はないと思う。
#時系列データ分析で単位根検定（ＡＤＦ検定）の結果、単位根あり、つまり,
#非定常過程であるという帰無仮説を棄却できなければ、暫定的に非定常過程と
#仮定して分析をすすめることになるが何となく落ち着きが悪い。実務的観点か
#ら分析者の主目的が予測モデルの診断にあるとすれば、原系列のままのlevel
#データと差分データを使ったモデルの２者を比べてみてout-of- sample テスト
#での予測精度よいほうのモデルを選択するのも一つの便法かもしれない。