参考：
http://www.jaceksuda.com/teach/pse2011/class9_slides.pdf
単位根検定→定常性の検定→RW性の検定
70年代後半まで：
線形トレンドでlog(GDP)をフィッティング、確率部分はここからのゆらぎ
70年代問題
log(GDP)が線形増加がスローダウン→予見できなくなる。
多項回帰などより複雑なトレンドが考えられるようになる
→低周波数領域では有効
ＧＮＰのショックはランダムウォークに似ている。
→1982年 Nelson and Plosserが単位根検定
→単位根の存在を棄却できず。
実際には、単位根を検定する
ARMAを利用
とりあえず、AR(1)の例：
wikiより，
http://ja.wikipedia.org/wiki/%E8%87%AA%E5%B7%B1%E5%9B%9E%E5%B8%B0%E7%A7%BB%E5%8B%95%E5%B9%B3%E5%9D%87%E3%83%A2%E3%83%87%E3%83%AB
x(t+1)=φx(t)+c+ε
→t→∞で
x(t)=c/(1-φ)+Σφ^kε[t-k]
φ^n→0 だと定常へ
ARMAは力を過去から推定していることに相当。現在力のノイズ消し。
ラグ演算子
ε(t)=(1-Σφ_iL_i)X(t)=φx(t)
φ=1-Σφ_iL_i （自分-過去の影響の演算子）
AR：φX(t)=ε(t)
MA: X(t)=(1+Σ_{i=1 to q}θ_iL_i=θε_(t)
ここでθ=Σ_{i=1 to q}θ_iL_i←過去のショックと今のショックの合計（外力演算子）
したがってARMAは
φX(t)=θε(t)
→線形代数的に、
φ(z)=0 の単位根が１だと定常にならない。
→ということで単位根＝１を検定
DF検定：場合わけがある
(1)単純なランダムウォークの場合
x(t+1)=φx(t)+ε
φを線形回帰で推定Ａ
t_{φ}=(A-φ)/Se(A)を検定
もしφが0.9とか0.8とかなら、tはt分布
しかし、φ=1のときは解析的に単純にかけない分布に収束
→DF分布。DF分布はウィナー過程の積分的なもの。
→DF分布で検定。
(2)Constant+AR(1)
(y(t+1)-μ)=φ(y(t-1)-μ)+ε(t)
t^{μ}_{φ}というμに依存する検定量をつかう。
これはμ=1のときDF^{μ}分布に従う。
(３)ドリフト+定数+ランダムウォーク
y(t)-c-β(t)=φ(y(t-1)-c-β(t-1))+ε(t)
→t^{β}_{φ=1}統計量はDF^{β}分布に従う。
β=0のときは、DF^{μ}の分布のが強力な検定。
ADFtest(1984) DFtestの拡張版ARMAに対応
系列相関がある場合
Cov(u_i,u_j)≠0 ということ
Phillips-Perron Unit-Root Test：どんな系列相関に対応
検定量
Z(t)=(σ^2/λ^2)t_{ρ=0}-1/2*1(T・SE(ρ)/σ^2)
t_{ρ=0}=ρ/SE(ρ)
σ2^とλ^2の推定
σ^2=T^{-1}ΣE(u^2)
λ^2=ΣE(T^-1S_T^2)
S_T=Σu_t
u_tは誤差項、
Under the null hypothesis （帰無仮説 ρ=0 ←ん？、帰無仮説はρ=1?）では、tはＡＤＦ分布に従う。
ほかのλの推定
λ^2=γ0+2Σ(1-j/(m+1))γ_j
γ_0=TΣu(t)^2
γ_j=1/TΣ_{t=j+1? to T}u_tu_{t-j}
http://faculty.washington.edu/ezivot/econ584/notes/unitroot.pdf
http://economia.unipv.it/pagp/pagine_personali/erossi/rossi_unit_roots_PhD.pdf
ssqrtl（λ^2）の計算
if (lshort)
l <- trunc(4 * (n/100)^0.25)
else l <- trunc(12 * (n/100)^0.25)
ssqrtl <- .C("R_pp_sum", as.vector(u, mode = "double"), as.integer(n),
as.integer(l), trm = as.double(ssqru), PACKAGE = "stats")
ssqrtl <- ssqrtl$trm
lshortはデフォルトがTらしい
.CはCの関数を実行。
R_pp_sum はおそらくこれ
http://fossies.org/dox/R-2.15.2/PPsum_8c_source.html
で、lが謎杉
実質やっていることは、nは回帰のデータ数、lは上で決まる定数
gamma2<-rep(0,l);for(j in 1:l){for(t in (j+1):n){gamma2[j]<-gamma2[j]+u[t-j]*u[t];}};gamma2<-gamma2/n
gamma2<-gamma0+2*sum([1-(1:l)/(1+l)]*gamma2)
gamma2
この部分 tableでいくつかのパタンで近似
tableは25,50,100,250,500,10^5の値が入ってる。
例えば、一つ目だと4.38,4.15,4.04,3.99,3.96
各パタンについて、それぞれnのときの値を計算．
approx(x,y,n,rule=2)
xとyの値を内挿して, nの値を求める、かつ, rule=2とする。
ルール２：外装状態のときの処理：最大値を使う (n=10^8 のときは3.96になる）
ルール１：計算しないNAにする
table <- cbind(c(4.38, 4.15, 4.04, 3.99, 3.98, 3.96), c(3.95,
3.8, 3.73, 3.69, 3.68, 3.66), c(3.6, 3.5, 3.45, 3.43,
3.42, 3.41), c(3.24, 3.18, 3.15, 3.13, 3.13, 3.12), c(1.14,
1.19, 1.22, 1.23, 1.24, 1.25), c(0.8, 0.87, 0.9, 0.92,
0.93, 0.94), c(0.5, 0.58, 0.62, 0.64, 0.65, 0.66), c(0.15,
0.24, 0.28, 0.31, 0.32, 0.33))
table <- -table
tablen <- dim(table)[2L]
tableT <- c(25, 50, 100, 250, 500, 1e+05)
tablep <- c(0.01, 0.025, 0.05, 0.1, 0.9, 0.95, 0.975, 0.99)
tableipl <- numeric(tablen)
for (i in (1L:tablen)) tableipl[i] <- approx(tableT, table[,
i], n, rule = 2)$y
*1:λ^2-σ^2)/(λ^2
