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

#ディッキー(Dickey)とフラー(Fuller)はランダムウォーク過程では最小
#2乗推定量にバイアスが生じることを示した。
#単位根検定(unit root test)(ADF検定)ディッキーフラー(Dickey-Fuller)
#検定の考え方と背景について簡単なシミュレーションとグラフで検討する。
#まずはじめに、ランダムウォーク過程の分散について見てみる。

#
ランダムウォーク過程の分散
#ランダムウォーク過程に従う時系列データの特徴をグラフで
#見てみる。u(t)は期待値0、分散σ² のホワイトノイズとする。
#ホワイトノイズは系列相関がゼロを前提にするが、i.i.d(independently
# and identically distributed)のようにu(t)が独立であることまで
# 要求しないという意味で、i.i.dより弱い条件となっている。(相関が
#ゼロであったとして必ずしも独立とは限らないという意味で)
#ランダムウォーク過程を y(t)=y(t-1)+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 となることが、どの時系列分析のテキスト
#にも書かれている。ランダムウォーク過程に従う時系列はその分散が一定
#でなく時間の経過と共に増加していく特徴がある。
#この現象を正規分布(平均値ゼロ、標準偏差 1)に従うホワイトノイズ
#によるシミュレーションで表すことを試みる。
#手順としては100期間(例えば時間であ100時間)の時系列を
#2000本作成してその時間の経過後との平均値と分散の変化を視覚化して
#みる。
#

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

Nt<-100  # サンプル数(期間数)

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

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

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


#2000本のランダムウォーク系列を作成する
rwalk1<-apply(rwalk,2,cumsum)
#2000本のランダムウォーク系列のグラフ
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乗推定量のバイアスについて
#このように分散が時間の経過と共に増加する時系列データに
#最小二乗法を適用してパラメータを推定した場合にはバイアスが生じる
#ことが知られている。回帰分析などを行った場合には時間の経過と共に
#分散が変化するのでt分布の形も変わってしまう。
#ディッキー(Dickey)とフラー(Fuller)はランダムウォーク過程では最小2乗
#推定量にバイアスが生じることを示した。
#真のモデルがy(t)=ρ・y(t-1)とした場合に
#両辺からy(t-1)を差し引いて
#y(t)-y(t-1)=△y(t),β=ρ-1 とすれば
# △y(t)=β・y(t) で表現できる。この式について
#最小二乗法を適用してβを推定するとρが1であれば、つまり単位根が存在
#すればβはゼロとなることが予想される。
#
#そこで3種類のランダムウォーク
#(定数項(ドリフト α0)が無い場合
#△y(t)=β・y(t)+ u(t)
#ドリフト α0 がある場合
#△y(t)=α0 + β・y(t)+ u(t)
#、トレンド項tとドリフトα0の両者がある場合
#△y(t)=α0 + β・y(t)+ α1・t+u(t)
#の三種類のケースについて最小2乗法を使って回帰係数とt値を推定してみる。

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

bvalnoconst<-bvaldrift<-bvaltrend<-rep(0,nt) 
tvalnoconst<-tvaldrift<-tvaltrend<-rep(0,nt) 
#下記でyはランダムウォーク系列y(t)を示し
#yyはy(t-1) つまり1期遅れの系列を示し
#dyは△y(t)を示している。
#あとは回帰分析の結果(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="ドリフトなし βの分布")

ドリフトなし βの分布  単位根過程


#真のモデルがランダムウォーク過程に従っていれば、ρ=1つまり単位根が

# 存在するのでβ=ρ-1=0 となる筈であるがβの推定値はゼロより小さい
#マイナス側に偏っていることが分かる。ドリフト付き、トレンド付きになると
#マイナス側への偏りは、さらに大きくなる。仮にβ=-0.1となれば
#ρ=0.9であり、単位根が真に存在するときでも最小2乗推定値は0.9
#と計算する可能性を意味している。

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

ドリフトつき  βの分布  単位根過程

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

#3種類もモデルについてのt値の分布をグラフにすると左右対照のt分布

#とは異なっておりマイナスである側に偏っていることがわかる。
#単位根の検定をするときには通常のt分布表を使うことは適当でなく
#マイナス側への偏りを考慮する必要があることがわかる。ざっくり言えば
#より大きいマイナスのt値が出てくれば帰無仮説の単位根あり(非定常)
#を棄却できることがわかる。


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

ドリフトなし  t値の分布  単位根検定 DF分布


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

ドリフト付き  t値の分布  単位根検定 DF分布



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

ドリフト トレンド付き  t値の分布  単位根検定 DF分布


print("ドリフトなし  サンプル数 100")
## [1] "ドリフトなし  サンプル数 100"
quantile(tvalnoconst,c(0.01,0.05,0.1))
##      1%      5%     10%

## -2.5517 -1.9384 -1.6248
print("ドリフト付き  サンプル数 100")
## [1] "ドリフト付き  サンプル数 100"
quantile(tvaldrift,c(0.01,0.05,0.1))
##      1%      5%     10%

## -3.6705 -2.9017 -2.6580
print("ドリフト、トレンド 付き  サンプル数 100")
## [1] "ドリフト、トレンド 付き  サンプル数 100"
quantile(tvaltrend,c(0.01,0.05,0.1))
##      1%      5%     10%

## -4.0456 -3.4542 -3.1615
#βの推定値の分布はゼロより少し偏小さい数値に偏って分布している。

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

 

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

 

金融工学の初等的計算例 目次