読者です 読者をやめる 読者になる 読者になる

マツリさんの日記

androidプログラミング初心者の奮闘日記です。たまに統計学もしてます。

時系列データ解析 06 ~ 自己回帰過程

www.asakura.co.jp

 沖本先生の本の続きです。MA過程の次はAR過程です。

 自己回帰(AR)過程は、過程が自己の過去に回帰されたモデルで表現されるものです。1次のAR過程、AR(1)は、

 \displaystyle y_{t} = c + \phi_{1} y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{2})

で定義され、 y_{t} \sim AR(1) と表記されます。

 AR過程は代表的な弱定常過程である

 \displaystyle y_{t} = \mu + \epsilon_{t}

に、 \phi_{1} y_{t-1} が追加された形をしていて、そのため y_{t} y_{t-1} が相関をもつことになります。

 AR過程の確率的変動は、MA過程と同様、撹乱項のホワイトノイズ  \epsilon_{t} によって決まります。 \epsilon が決まり、続いて y が決まっていきます。

 AR過程では初期値の決定が問題になります。 y が分布が定まっている場合は、その分布に従う確率変数とし、 y の分布が定まっていない場合は、ある定数とすることが一般的です。ただ、定常過程では、初期値の影響は時間とともになくなるので、大きな問題にはなりません。

 沖本先生の本のとおり、具体例を示します。次のようなAR(1)を考えます。

 y_{t} = 1 + 0.5 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N.(1)

 初期値を  y_{0} = 0 として、 \epsilon_{1} = - 2.1 \epsilon_{2} = 0.9 とします。このときに順番に解いていくと、

 y_{1} = 1 + 0.5 y_{0} \epsilon_{1} = 1 + 0.5 \dot - 2.1 = - 1.1
 y_{2} = 1 + 0.5 y_{1} \epsilon_{2} = 1 + 0.5 \dot (-1.1) + 0.9 = 1.35

となります。

 ここで、MA(1) 過程と同様に、パラメータを調整したAR(1)過程をグラフ化してみます。RのコマンドについてはMA(1)で多少説明していますので、必要な部分のみ解説します。

 以下の6つのパラメータのAR(1)過程をグラフ化します。

 \displaystyle \begin{align}
&(a) y_{t} = 2 + 0.8 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{1}) \\ 
&(b) y_{t} = - 2 + 0.3 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{0.5}) \\
&(c) y_{t} = 0 - 0.3 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{2}) \\
&(d) y_{t} = - 2 - 0.8 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{1}) \\
&(e) y_{t} = 2 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{1}) \\
&(f) y_{t} = 2 + 1.1 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{1}) \end{align}

 Rのコードです。arima.sim関数の第2引数ですが、今回はAR過程なので、order=c(1, 0, 0) となっています。ちなみに、MA過程の時はorder=c(0, 0, 1) でした。次にAR過程の係数はar=c(0.8)として指定します。

set.seed(1)
d1 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(0.8)), sd=1) +2
plot(d1, main="(a) c=2, Φ=0.8, σ=1")
set.seed(1)
d2 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(0.3)), sd=0.5) - 2
plot(d2, main="(b) c=-2, Φ=0.3, σ=0.5")
set.seed(1)
d3 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(-0.3)), sd=2)
plot(d3, main="(c) c=0, Φ-=0.3, σ=2")
set.seed(1)
d4 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(-0.8)), sd=1) - 2
plot(d4, main="(d) c=-2, Φ=-0.8, σ=1")
set.seed(1)
d5 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(1)), sd=1) + 2
plot(d5, main="(e) c=2, Φ=1, σ=1")
set.seed(1)
d6 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(1.1)), sd=1) + 2
plot(d6, main="(f) c=2, Φ=1.1, σ=1")

 グラフ化するとこのようになります。


 まず、(e)と(f)ですが、係数が 1以上の場合、AR過程は定常ではなくなるため、エラーで出力されません。そのため、グラフも出力されません。

 AR過程はMA過程と違い、パラメータの値によって定常か非定常かが変わります。先ほど、係数が 1以上と書きましたが、正確には係数 |\phi_{1}| \ge 1の時に非定常となります。
 ここではAR(1)過程が定常という場合のみを解説していきます。

 AR(1)過程の期待値についてです。グラフでは概ね期待値と定数項が一致しているようにも見えますが、実際は期待値と定数項は一致しません。(a)あたりでは定数項は 2ですが、やや下の方に過程の期待値があるようです。AR(1)過程の期待値をみています。

 \displaystyle y_{t} = c + \phi_{1} y_{t-1} + \epsilon_{t}
 \displaystyle E(y_{t}) = E(c + \phi_{1} y_{t-1} + \epsilon_{t}) = c + \phi_{1} E(y_{t-1})

ここで、ホワイトノイズの性質により E(\epsilon_{t}) = 0 です。また、 yは定常過程ですから、

 \displaystyle E(y_{t}) = E(y_{t-1}) = \mu

です。よって、AR(1)過程の期待値

 \displaystyle \mu = c + \phi_{1} \mu

 \displaystyle \mu = \frac{c}{(1 - \phi_{1})}

が得られます。

 次に分散ですが、MA過程と同様に過程の分散は撹乱項の分散 \sigma^{2}と異なります。MA(1)過程で使用した(a)のグラフと今回の(a)のグラフを比べてみます。ちなみに左側はMA(1)過程、右側がAR(1)過程です。


 撹乱項の分散は 1ですが、過程の変動幅は違います。AR(1)過程の方がやや大きいです。つまり、AR(1)過程の分散は撹乱項のそれより大きくなります。では分散を見ていきます。

 \displaystyle y_{t} = c + \phi_{1} y_{t-1} + \epsilon_{t}
 \displaystyle \begin{align} Var(y_{t}) &= Var(c + \phi_{1} y_{t-1} + \epsilon_{t}) \\
&= \phi_{1}^{2} y_{t-1}^{2} + 2 \phi_{1} y_{t-1} (c + \epsilon_{t} - \mu) + (c + \epsilon_{t} - \mu)^{2} \\
&= \phi_{1}^{2} Var(y_{t-1}) + 2 \phi_{1} Cov(y_{t-1} , \epsilon_{t}) + Var(\epsilon_{t}) \\
&= \phi_{1}^{2} Var(y_{t-1}) + \sigma^{2} \end{align}

 ここで、 y_{t}が定常であれば、 \displaystyle \gamma_{0} = Var(y_{t}) = Var(y_{t-1})ですので、AR(1)過程の分散は、
 \displaystyle \gamma_{0} = \frac{\sigma^{2}}{(1 - \phi_{1}^{2})}
となります。

 最後に自己相関をみてみましょう。先ほどのMA(1)過程のグラフ(a)とAR(1)過程のグラフを比べると、係数( \theta_{1} \phi_{1} はともに 0.8)は同じでもAR(1)過程のグラフがより滑らかであることがわかります。

 これは、この場合のように係数が同じとき、1次自己相関の絶対値はAR(1)過程の方が大きくなることと、AR(1)過程は2次以降の相関係数もすべて正になることの2点によるからです。

 AR過程のk次の自己共分散はこのようになります。

 \displaystyle \begin{align} \gamma_{k} &= Cov(y_{t} , y_{t-k}) \\
&= Cov(\phi_{1} y_{t-1} + \epsilon_{t} , y_{t-k}) \\
&= Cov(\phi_{1} y_{t-k} , y_{t-k}) + Cov(\epsilon_{t} , y_{t-k}) \\
&= \phi_{1} \gamma_{k-1} \end{align}

 ここで、両辺を \gamma_{0} で割ると、

 \displaystyle \gamma_{k} = \phi_{1} \rho_{k-1}

となり、これはユール・ウォーカー方程式と呼ばれます。ユール・ウォーカー方程式は、AR過程の自己相関が、 y_{t} が従うAR過程と同一の係数 \phi_{1} をもつ差分方程式に従うということを表しています。

 ユール・ウォーカー方程式と \rho_{0} = 1を用いて、AR過程の自己相関は順番に求めることができます。また、ユール・ウォーカー方程式により、AR(1)過程の一般的な自己相関は、

 \displaystyle \rho_{k} = \phi_{1}^{2} \rho_{k-2} = \phi_{1}^{3} \rho_{k-3} = \cdots = \phi_{1}^{k} \rho_{0} = \phi_{1}^{k}

となります。また、[ | \phi_{1} | < 1 より、AR(1)過程の自己相関の絶対値は指数的に減衰していきます。コレログラムを描いてみます。

set.seed(1)
d1 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(0.8)), sd=1)
acf(d1, main="(a)  θ1=0.8, θ2=0")
set.seed(1)
d2 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(-0.8)), sd=1)
acf(d2, main="(a)  θ1=0.8, θ2=0")

 コレログラムを描くacf()関数とarima.sim()関数は以前に説明していますのでここでは割愛します。AR(1)過程の自己相関の絶対値は指数的に減衰していくことが、(a)と(b)からよくわかると思います。AR(1)過程のコレログラムは \phi_{1}の符号によって、減衰していくことが確認できます。

 AR(1)過程を一般化したp次のAR過程AR(p)はこのように記述されます。

 \displaystyle y_{t} = c + \phi_{1} y_{t-1} + \phi_{2} y_{t-2} + \cdots + \phi_{p} y_{t-p} + \epsilon_{t}

 AR(p)過程はy_{t} を定数と自身のp期過去の値に回帰したモデルになります。AR(1)同様にAR(q)もパラメータにより定常となるか非定常となるか変わってきます。定常AR(p)過程の性質は、

 \displaystyle \mu = E(y_{t}) = \frac{c}{1 - \phi_{1} - \phi_{2} - \cdots - \phi_{p}}
 \displaystyle \gamma_{0} = Var(y_{t}) = \frac{\sigma^{2}}{1 - \phi_{1} \rho_{1} - \phi_{2} \rho_{2} - \cdots - \phi_{p} \rho_{p}}
 \displaystyle \gamma_{k} = \phi_{1} \gamma_{k-1} + \phi_{2} \gamma_{k-2} + \cdots + \phi_{p} \gamma_{k-p}, k \ge 1
 \displaystyle \rho_{k} = \phi_{1} \rho_{k-1} + \phi_{2} \rho_{k-2} + \cdots + \phi_{p} \rho_{k-p}, k \ge 1

であり、特に \gamma_{k}はユール・ウォーカー方程式と呼ばれます。そして、自己相関が指数的に減衰していきます。

 ここで、AR(2)過程のユール・ウォーカー方程式は、このようになります。

 \displaystyle \rho_{k} = \phi_{1} \rho_{k-1} + \phi_{2} \rho_{k-2}

これに、 k=1を代入し、 \rho_{0} = 1 \rho_{1} = \rho_{-1}から

 \displaystyle \rho_{1} = \phi_{1} \rho_{1} + \phi_{2} \rho_{-1} = \phi_{1} + \phi_{2} \rho_{1}

が得られます。これを変形すると、

 \displaystyle \rho_{1} = \frac{\phi_{1}}{(1 - \phi_{2})}

となります。次に、 k=2を代入すると、

 \displaystyle \rho_{2} = \phi_{1} \rho_{1} + \phi_{2} \rho_{0} = \frac{\phi_{1}^{2}}{1 - \phi_{2}} + \phi_{2} = \frac{\phi_{1}^{2} + \phi_{2} - \phi_{2}^{2}}{1 - \phi_{2}}

となり、 \rho_{2} を求めることができます。AR(2)過程の場合は、特に \rho_{k}の一般的な

 \displaystyle 1 - \phi_{1} z - \phi_{2} z^{2} = 0

という2次方程式が異なる2つの解をもつとき、その解の逆数を \lambda_{1} \lambda_{2}とすれば

 \displaystyle \rho_{k} = \frac{(1 - \lambda_{2}^{2}) \lambda_{1}^{k+1} - (1 - \lambda_{1}^{2}) \lambda_{2}^{k+1}}{(\lambda_{1} - \lambda_{2})(1 + \lambda_{1} \lambda_{2})}

で自己相関係数が得られます。

 AR(p)の場合、 \rho_{1}から \rho_{k-1}までは、ユール・ウォーカー方程式と \rho_{0} = 1 \rho_{k} = \rho_{-k}から p-1次の連立方程式を解いて求めることができます。

 AR(2)過程のコレログラムを図示します。

set.seed(1)
d3 <- arima.sim(n=100, model=list(order=c(2, 0, 0), ar=c(0.5, 0.35)), sd=1)
acf(d3, main="(a)  θ1=0.5, θ2=0.35")
set.seed(1)
d4 <- arima.sim(n=100, model=list(order=c(2, 0, 0), ar=c(0.1, 0.5)), sd=1)
acf(d4, main="(a)  θ1=0.1, θ2=0.5")
set.seed(1)
d5 <- arima.sim(n=100, model=list(order=c(2, 0, 0), ar=c(0.5, -0.8)), sd=1)
acf(d5, main="(a)  θ1=0.5, θ2=-0.8")
set.seed(1)
d6 <- arima.sim(n=100, model=list(order=c(2, 0, 0), ar=c(0.9, -0.8)), sd=1)
acf(d6, main="(a)  θ1=0.9, θ2=-0.8")

 AR(1)同様に自己相関の絶対値は指数的に減衰しています。パラメータの組み合わせによって、多様な自己相関構造が記述できます。(e)と(f)をみると、AR(2)が循環的な自己相関構造を記述できることがわかると思います。先ほどの

 \displaystyle 1 - \phi_{1} z - \phi_{2} z^{2} = 0

が共役な複素数 a \pm b \mathrm{i}を解にもつときに、自己相関は循環的になります。そして、その周期は

 \frac{2 \pi}{\cos^{-1} \left[ a / \sqrt{a^{2} + b^{2}}\right]}

となります。これは、AR(p)過程にも拡張でき、

 \displaystyle 1 - \phi_{1} z - \cdots - \phi_{p} z^{p} = 0

が共役な複素数を解にもつときにも同様の循環成分をもちます。AR(p)の自己相関は、最大 \frac{p}{2}個の循環成分を持つことができるのです。

時系列データ解析 05 ~ 移動平均過程

www.asakura.co.jp


 沖元先生の本の続きです。ここでは、移動平均モデルについての話からです。

 まず、1次MA過程(以下、MA(1)過程)はこのようにモデル化されます。

 \displaystyle y_{t} = \mu + \epsilon_{t} + \theta_{1} \epsilon_{t - 1}, \epsilon \sim W.N. (\sigma^{2})

 移動平均過程は、ホワイトノイズの線形和になります。

 ホワイトノイズはすべての時点において期待値が  0 で、分散が一定の確率過程です。また、自己相関をもたないので、弱定常過程となります。

  y_{t} が、移動平均過程に従うことを、  y_{t} \sim MA(1) とあらわします。

 もっとも基本的な弱定常過程は

 \displaystyle y_{t} = \mu + \epsilon_{t}, \epsilon \sim W.N. (\sigma^{2} )

ですから、移動平均過程はこの過程に、  \theta_{1} \epsilon_{t - 1} という項が追加されたものになります。

 また、 t期のMA(1)過程と  t - 1 期のMA(1)を比較してみます。

 \displaystyle y_{t} = \mu + \epsilon_{t} + \theta_{1} \epsilon_{t - 1}
 \displaystyle y_{t - 1} = \mu + \epsilon_{t - 1} + \theta_{1} \epsilon_{t - 2}

 この二つのモデルは、  \epsilon_{t - 1} という項を共通にもちます。そのため、  y_{t} y_{t - 1} に相関が生じます。

 そして、MA(1)過程は、弱定常過程と比べて、 \theta_{1} というパラメータが多くなっています。このパラメータが、自己相関の強さを決定づけます。

 MA(1)過程の確率変動は、ホワイトノイズ  \epsilon_{t} によって決まります。まず、  \epsilon が決まり、続いて  y が決まっていきます。沖本先生の本では、具体例として、

 \displaystyle y_{t} = 1 + \epsilon_{t} + 0.5 \epsilon_{t - 1}, \epsilon_{t} \sim W.N.(1)

というMA(1)過程で、  \epsilon_{0} = 0.5 \epsilon_{1} = - 2.1 \epsilon_{2} = 0.9 とすれば、

 \displaystyle y_{1} = 1 + \epsilon_{1} + 0.5 \epsilon_{0} = 1 - 2.1 + 0.5 \cdot 0.5 = - 0.85

 \displaystyle y_{2} = 1 + \epsilon_{2} + 0.5 \epsilon_{1} = 1 + 0.9 + 0.5 \cdot ( - 2.1 ) = - 0.85

という順番で  y が決められることが示されています。

 ここで、  \epsilon を正規ホワイトノイズとして、それぞれのパラメータの組み合わせたMA(1)をグラフ化してみます。

 Rを使って、シミュレーションをしてみます。Rにはさまざまなツールが用意されていまして、ここではARIMAを例にMA(1)過程をみていきます。ちなみに、ARIMA過程はこのようなコードで記述します。ARIMA過程そのものについては次回以降に説明します。

d <- arima.sim(n=100, model=list(order=c(2, 0, 4), ar=c(0.1, 5), ma=c(-0.2, 0.3, 0.4, -0.5)), sd=1)

 arima.simの引数について説明します。n=100がシミュレーションの数です。今回の話の流れでは、時系列のデータを  T = 100 まで実験してみるということです。

 第2引数は、model=list()となっていますが、中身をみていきます。order=()の引数は、自己回帰過程、和分過程、移動平均過程の次数になります。自己回帰過程、和分過程は後述します。ar=()の引数は、係数で、ma=()も同様です。

 第3引数については、分散の大きさといいますか、上の例では分散も標準偏差(標準誤差)も  1です。

 それでは、シミュレーションしてみます。沖本先生の例にならって6パターンにします。

 \displaystyle \begin{align}
&(a) y_{t} = 0 + \epsilon_{t} + 0.8 \epsilon_{t - 1}, \epsilon \sim W.N. (\sigma^{1}) \\
&(b) y_{t} = 2 + \epsilon_{t} + 0.5 \epsilon_{t - 1}, \epsilon \sim W.N. (\sigma^{0.5}) \\
&(c) y_{t} = - 2 + \epsilon_{t} + 0.3 \epsilon_{t - 1}, \epsilon \sim W.N. (\sigma^{2}) \\
&(d) y_{t} = 0 + \epsilon_{t} - 0.3 \epsilon_{t - 1}, \epsilon \sim W.N. (\sigma^{1}) \\
&(e) y_{t} = 2 + \epsilon_{t} - 0.5 \epsilon_{t - 1}, \epsilon \sim W.N. (\sigma^{0.5}) \\
&(f) y_{t} = - 2 + \epsilon_{t} - 0.8 \epsilon_{t - 1}, \epsilon \sim W.N. (\sigma^{2}) \end{align}

 ついでにコードも記述します。

set.seed(1)
d1 <- arima.sim(n=100, model=list(order=c(0, 0, 1), ma=c(0.2)), sd=1)
plot(d1, main="(a) μ=0, θ=0.8, σ=1")
set.seed(1)
d2 <- arima.sim(n=100, model=list(order=c(0, 0, 1), ma=c(0.5)), sd=0.5) + 2
plot(d2, main="(b) μ=2, θ=0.5, σ=0.5")
set.seed(1)
d3 <- arima.sim(n=100, model=list(order=c(0, 0, 1), ma=c(0.3)), sd=2) - 2
plot(d3, main="(c) μ=-2, θ=0.3, σ=2")
set.seed(1)
d4 <- arima.sim(n=100, model=list(order=c(0, 0, 1), ma=c(-0.3)), sd=1)
plot(d4, main="(d) μ=0, θ=-0.3, σ=1")
set.seed(1)
d5 <- arima.sim(n=100, model=list(order=c(0, 0, 1), ma=c(-0.5)), sd=0.5) + 2
plot(d5, main="(e) μ=2, θ=-0.5, σ=0.5")
set.seed(1)
d6 <- arima.sim(n=100, model=list(order=c(0, 0, 1), ma=c(-0.8)), sd=2) - 2
plot(d6, main="(f) μ=-2, θ=-0.8, σ=2")

 ar=()については省略しています。plot()はグラフに表示するコマンドです。main=""はグラフのタイトルのようなものです。set.seet(1)は、乱数種を指定して乱数を発生させる関数です。すべて同じ乱数を発生させてシミュレーションしています。

 グラフ化するとこうなります。



 グラフから確認できることは、系列が  \mu の周りを変動していることです。つまり、MA(1)過程の期待値は  \mu になります。

 \begin{align} E(y_{t}) &= E(\mu + \epsilon_{t} + \theta_{1} \epsilon_{t - 1}) \\
&= E(\mu) + E(\epsilon_{t}) + E(\theta_{1} \epsilon_{t - 1}) \\
&= E(\mu) + E(\epsilon_{t}) + \theta_{1} E(\epsilon_{t - 1}) \\
&= \mu \because E(\epsilon_{t}) = 0, E(\epsilon_{t - 1}) = 0)\end{align}

 次にMA(1)過程の分散は、撹乱項  \epsilon_{t} の分散  \sigma^{2} に等しいと既に説明しました。しかし、グラフを見ると分かるとおり、撹乱項の分散が同じでも若干変動幅が異なることが分かると思います。

 つまり、MA(1)過程分散は、撹乱項の分散よりも大きくなるのです。

 \displaystyle y_{t} - \mu = \epsilon_{t} - \theta_{1} \epsilon_{t - 1}

 \begin{align} (y_{t} - \mu)^{2} &= (\epsilon_{t} - \theta_{1} \epsilon_{t - 1})^{2} \\
&= \epsilon_{t}^{2} - 2 \theta_{1} \epsilon_{t} \epsilon_{t - 1} + \theta_{1}^{2} \epsilon_{t - 1}^{2} \end{align}

 \begin{align}
\gamma_{0} &= E \left[ \epsilon_{t}^{2} \right] - 2 E \left[ \theta_{1} \epsilon_{t} \epsilon_{t - 1} \right] + \theta_{1}^{2} E \left[ \epsilon_{t - 1}^{2} \right] \\
&= (1 + \theta_{1}^{2} ) \sigma^{2} \end{align}

 したがって、MA(1)過程の分散は \theta_{1}^{2} \sigma^{2} だけ、撹乱項の分散よりも大きくなります。

 また、グラフの滑らかさは  \theta_{1} が大きくなるにしたがって、より滑らかになります。(a)が一番滑らかであり、(f)が一番ギザギザしているます。これは、MA(1)過程の自己相関が  0 でないためです。MA(1)過程をもう一度見てみましょう。

 \displaystyle y_{t} = \mu + \epsilon_{t} + \theta_{1} \epsilon_{t - 1}, \epsilon
 \displaystyle y_{t - 1} = \mu + \epsilon_{t - 1} + \theta_{1} \epsilon_{t - 2}

  \theta_{1} が 正の場合、 y_{t} \theta_{1} \epsilon_{t-1} y_{t-1} \epsilon_{t-1} が同じ正の符号になるため、 y_{t} y_{t-1} が同じ方向に動く傾向をもちます。

 その結果、1次の正の自己相関が生じて、その自己相関は  \theta_{1} の値が  1 に近づくほど強くなります。

 逆に、 \theta_{1} が負の場合は、  y_{t} y_{t-1} は逆に動く傾向があり、1次の負の自己相関が生まれます。そのため、グラフはよりギザギザになります。

 次に、MA(1)過程の自己共分散を求めます。

 \displaystyle \begin{align} \gamma_{1} &= Cov( y_{t}, y_{t-1}) \\
&= Cov(\mu + \epsilon_{t} + \theta_{1} \epsilon_{t-1},  \mu + \epsilon_{t-1} + \theta_{1} \epsilon_{t-2} \\
&= Cov(\epsilon_{t},  \epsilon_{t-1}) + Cov(\epsilon_{t},  \theta_{1} \epsilon_{t-2} ) + Cov(\theta_{1} \epsilon_{t-1} , \epsilon_{t-1}) + Cov(\theta_{1} \epsilon_{t-1}, \theta_{t-1} \epsilon_{t-2}) \\
&= \theta_{1} Cov(\epsilon_{t-1}, \epsilon_{t-1}) \\
&= \theta_{1} \sigma_{2} \end{align}

 これは、ホワイトノイズ \epsilon k 次の自己共分散が  0 であり、分散が  0 であることから導けます。

 よって、MA(1)過程の1次自己相関を求めると、

 \displaystyle \rho_{1} = \frac{\gamma_{1}}{\gamma_{0}} = \frac{\theta_{1}}{1 + \theta_{1}^{2}}

となります。MA(1)過程の自己相関が  \mp 1 のとき、 \frac{1}{2} になります。したがいまして、1次の自己相関の絶対値が \frac{1}{2} より大きいMA(1)過程はモデル化できません。

 MA(1)過程の2次以降の自己共分散は、

 \displaystyle \begin{align} \gamma_{k} &= Cov(y_{t},  y_{t-k}) \\
&= Cov(\mu + \epsilon_{t} + \theta_{1} \epsilon_{t-1}, \mu + \epsilon_{t-k} + \theta_{1} \epsilon_{t-k-1}) \\
&=0 \end{align}

となります。つまり、MA(1)過程の2次以降の自己相関は  0 ということです。そのため、MA(1)過程は1次の自己相関をモデル化できますが、2次以降の自己相関を記述することはできません。

 以上のとおり、期待値も分散も  t によらず、定常性の性質を満たすことがわかりました。

 最後に、MA(1)過程のコレログラムを確認してみます。まずRのコードを記載しておきます。

set.seed(1)
d1 <- arima.sim(n=100, model=list(order=c(0, 0, 2), ma=c(0.2, 0)), sd=1)
acf(d1, main="(a)  θ1=0.8, θ2=0")
set.seed(1)
d2 <- arima.sim(n=100, model=list(order=c(0, 0, 2), ma=c(-0.8, 0)), sd=1)
acf(d2, main="(b) θ1=-0.8, θ2=0")
set.seed(1)
d3 <- arima.sim(n=100, model=list(order=c(0, 0, 2), ma=c(0.8, 0.3)), sd=1)
acf(d3, main="(c) θ1=0.8, θ2=0.3")
set.seed(1)
d4 <- arima.sim(n=100, model=list(order=c(0, 0, 2), ma=c(0.3, 0.8)), sd=1)
acf(d4, main="(d) θ1=0.3, θ2=0.8")
set.seed(1)
d5 <- arima.sim(n=100, model=list(order=c(0, 0, 2), ma=c(0.8, -0.3)), sd=1)
acf(d5, main="(e) θ1=0.8, θ2=-0.3")
set.seed(1)
d6 <- arima.sim(n=100, model=list(order=c(0, 0, 2), ma=c(-0.8, -0.3)), sd=1)
acf(d6, main="(f)  θ1=-0.8, θ2=-0.3")

 acf()は自己相関関数、つまりコレログラムを表しています。先ほどのarima.sim()の中身ですが、今回はMA(1)過程のオーダーを2にしていますので、ma=c()中に指定する係数も2つになっています。


 まず、それぞれの係数 [tex \theta] を見るとコレログラムが係数に応じて変化しているのがわかると思います。(a)、(b)については、1次のみ有意な自己相関があり、2次以降は有意な自己相関はありません。

 一般化した q移動平均過程は、

 y_{t} = \mu + \epsilon_{t} + \theta_{1} \epsilon_{t-1} + \theta_{2} \epsilon_{t-2} + \cdots + \theta_{q} \epsilon_{t-q}, \epsilon_{t} \sim W.N.(\sigma^{2})

と表現でき、MA(q)過程と書きます。MA(q)過程は、現在とq期間の過去のホワイトノイズの線形和に定数 \mu を加えたものになります。

 MA(q)過程の性質は、

 \displaystyle E(y_{t}) = \mu

 \displaystyle \gamma_{0} = Var(y_{t}) = (1 + \theta_{1}^{2} + \theta_{2}^{2} + \cdots + \theta_{q}^{2}) \sigma^{2}

 \displaystyle \gamma_{k} = \begin{cases}{} (\theta_{k} + \theta_{1} \theta_{k+1} + \cdots + \theta_{q-k} \theta_{q}) \sigma^{2}, &1 \leq k \leq q \\ 0, & k \geq q + 1 \end{cases}

 \displaystyle \rho_{k} = \begin{cases}{} \frac{\theta_{k} + \theta_{1} \theta_{k+1} + \cdots + \theta_{q-k} \theta_{q}}{1 + \theta_{1}^{2} + \theta_{2}^{2} + \cdots + \theta_{q}^{2}}, & 1 \leq k \leq q \\ 0, & k \geq q + 1 \end{cases}

であり、重要なことはMA過程は定常過程であり、MA(q)過程のq+1次以降の自己相関は  0 ということです。したがいまして、q次の自己相関構造をモデル化するには、q個のパラメータが必要になります。長期間の自己相関のモデル化はあまりに多くのパラメータが必要になるため、望ましいことではありません。

 次回は自己相関過程の話に移ります。

PhoneStateListenerの作成 03

 前回PhoneStateListenerをつかって、着信状態に応じてToast表示する方法を実装しました。

 今回はMainActivity.javaは省略して、PhoneReceiver.javaをこのように実装しなおしました。

 

 PhoneReceiver.java

import android.content.Context;
import android.database.Cursor;
import android.provider.ContactsContract;
import android.telephony.PhoneNumberUtils;
import android.telephony.PhoneStateListener;
import android.telephony.TelephonyManager;
import android.widget.Toast;

public class PhoneReceiver extends PhoneStateListener {
private Context context;
PhoneReceiver(Context context) {
this.context = context;
}

// 通話状態の変化に応じて表示を変更する
@Override
public void onCallStateChanged(int state, String incomingNumber) {
super.onCallStateChanged(state, incomingNumber);

switch (state) {
// 着信時の処理内容
case TelephonyManager.CALL_STATE_RINGING:

Cursor addressTable = context.getContentResolver().query(ContactsContract.CommonDataKinds.Phone.CONTENT_URI, null, null, null, null);
if (addressTable != null) {
while (addressTable.moveToNext()) {
String phoneNumber = addressTable.getString(addressTable.getColumnIndex(ContactsContract.CommonDataKinds.Phone.NUMBER));
if (PhoneNumberUtils.compare(phoneNumber, incomingNumber)) {
Toast.makeText(context, "OK!" + incomingNumber, Toast.LENGTH_LONG).show();
}
}
addressTable.close();
}
break;

// 通話時の処理内容
case TelephonyManager.CALL_STATE_OFFHOOK:
Toast.makeText(context, "通話中!!" + incomingNumber, Toast.LENGTH_LONG).show();
break;
}
}
}

 実は更新をしていなかった間に、色々変更に変更を重ねて、このかたちになりました。

 当初、ContentResolverを用いて、データを抽出する作業をしていましたが、問題はデータの保存でして。

 SQLite、realmで保存方法を試行錯誤しましたが、どうしても納得できるかたちにできず。

 最終的に、Stack Overflowを参考にすると、PhoneNumberUtilsという便利なAPIを見つけることができました。

 日本語の解説がほぼなかったので、見つけるのに難儀しました。

時系列データ解析 04 ~ 自己相関の検定

www.asakura.co.jp


 沖本先生の本の続きです。今回は自己相関の検定です。
 そもそも時系列データに自己相関構造がなければ、モデルを構築して、そのモデルから予測を行うことも難しくなります。
 そのため、自己相関の有無を調べることは、大変重要なことになります。

 定常過程での基本統計量は、時点に依存しないため、標本統計量からすぐに計算できます。

 期待値、自己共分散、自己相関係数は、それぞれ次のようになります。

 \displaystyle \overline{y} = \frac{1}{T} \sum^{T}_{t = 1} y_{t}

 \displaystyle \hat{\gamma}_{k} = \frac{1}{T} \sum^{T}_{t = k + 1} (y_{t} - \overline{y }) (y_{t-k} - \overline{y} ), k = 0, 1, 2, \cdots

 \displaystyle \hat{\rho}_{k} = \frac{\hat{\gamma}_{k}}{\hat{\gamma}_{0}}, k = 1, 2, 3, \cdots

 それぞれ、標本平均、標本自己共分散、標本自己相関係数と呼ばれます。

 この標本自己相関係数  \hat{\rho}_{k} を使って、仮説検定を行います。
 帰無仮説は  H_{0} : \rho_{k} = 0、つまり自己相関がないということです。
 対立仮説は、  H_{1} : \rho \neq 0で、自己相関が  0 ではないということです。

 特に、  y_{t} iid 系列の場合、標本自己相関係数  \hat{\rho}_{k} が、漸近的に平均  0 、分散  \frac{1}{T}正規分布に従うそうです。これについて、統計量を求めることになります。

 複数の  k に関して、自己相関係数  \rho_{k} 0 の場合を検定する場合は、帰無仮説は  H_{0} : \rho_{1} = \rho_{2} = \cdots = \rho_{m} = 0 となります。対立仮説は、  H_{1} : 少なくとも一つの  k \in \left[ 1, m \right]において  \rho_{k} = 0 となります。

 この検定は、かばん検定(portmanteau test)と呼ばれ、Ljung Box検定が有名です。
 Ljung Box検定の統計量は、

 \displaystyle Q(m) = T(T + 2) \sum^{m}_{k = 1} \frac{\hat{\rho}^{2}_{k}}{T - k} \sim \chi^{2} (m)

が一定の条件のもとで成立します。この検定では、小さい  m を選択すると高次の自己相関を見逃したり、大きい  m を選択すると検定の検出力が小さくなるなどの問題があります。
  m を選択する際の目安として、  m \approx log(T) がありますが、実務的には複数の  m に関して検定を行っているようです。

 今更ですが、正規分布確率密度関数

 \displaystyle f(x) = \frac{1}{\sqrt{2 \pi} \sigma} exp \left\{ - \frac{(x - \mu)^{2}}{2 \sigma^{2}} \right\}

と表現されます。
 これで第1章は終了です。

PhoneStateListenerの作成 02

 毎度、小出しにしていますが、PhoneStateListenerの続きです。

 前回の実装では、バックグラウンド時にはToastが表示されないという点で未完成でした。

 そこで、今回は他のサイトのBroadcastReceiver、IntentFilterを参考にして、バックグランドでも通話状態を捕捉するように実装しなおしました。

MainActivity.java

import android.content.Context;
import android.support.v7.app.AppCompatActivity;
import android.os.Bundle;
import android.telephony.PhoneStateListener;
import android.telephony.TelephonyManager;

public class MainActivity extends AppCompatActivity {

// 各フィールドの設定
PhoneReceiver phoneStateListener;
TelephonyManager manager;

@Override
protected void onCreate(Bundle savedInstanceState) {
super.onCreate(savedInstanceState);
setContentView(R.layout.activity_main);

// PhoneReceiverインスタンスの生成
phoneStateListener = new PhoneReceiver(this);
// TelephonyManagerインスタンスの生成(Context.TELEPHONY_SERVICEを指定)
manager = ((TelephonyManager) getSystemService(Context.TELEPHONY_SERVICE));
}

@Override
protected void onResume() {
super.onResume();
manager.listen(phoneStateListener, PhoneStateListener.LISTEN_CALL_STATE);
}
}

PhoneReceiver.java

import android.content.Context;
import android.telephony.PhoneStateListener;
import android.telephony.TelephonyManager;
import android.widget.Toast;

public class PhoneReceiver extends PhoneStateListener {
Context context;
public PhoneReceiver(Context context) {
this.context = context;
}

// 通話状態の変化に応じて表示を変更する
@Override
public void onCallStateChanged(int state, String incomingNumber) {
super.onCallStateChanged(state, incomingNumber);

switch (state) {
// 着信時の処理内容
case TelephonyManager.CALL_STATE_RINGING:
Toast.makeText(context, "着信中!!" + incomingNumber, Toast.LENGTH_LONG).show();
break;
// 通話時の処理内容
case TelephonyManager.CALL_STATE_OFFHOOK:
Toast.makeText(context, "通話中!!" + incomingNumber, Toast.LENGTH_LONG).show();
break;
}
}
}

PhoneStateReceiver.java

import android.content.BroadcastReceiver;
import android.content.Context;
import android.content.Intent;
import android.telephony.PhoneStateListener;
import android.telephony.TelephonyManager;

public class PhoneStateReceiver extends BroadcastReceiver {
// 各フィールドの定義
TelephonyManager manager;
PhoneReceiver phoneStateListener;
static boolean listener = false;

// intent情報を処理する
@Override
public void onReceive(Context context, Intent intent) {
// PhoneReceiverインスタンスの生成
phoneStateListener = new PhoneReceiver(context);
// TelephonyManagerインスタンスの生成
manager =((TelephonyManager)context.getSystemService(Context.TELEPHONY_SERVICE));

if(!listener) {
manager.listen(phoneStateListener, PhoneStateListener.LISTEN_CALL_STATE);
listener = true;
}
}
}

AndroidManifest.xml

<?xml version="1.0" encoding="utf-8"?>
<manifest xmlns:android="http://schemas.android.com/apk/res/android"
package="com.example.ma2ri.telephonecontrol">

<!-- 電話番号、通話の状態などの端末情報を取得するためのパーミッション -->
<uses-permission android:name="android.permission.READ_PHONE_STATE" />

<application
android:allowBackup="true"
android:icon="@mipmap/ic_launcher"
android:label="@string/app_name"
android:supportsRtl="true"
android:theme="@style/AppTheme">
<activity android:name="com.example.ma2ri.telephonecontrol.MainActivity">
<intent-filter>
<action android:name="android.intent.action.MAIN" />

<category android:name="android.intent.category.LAUNCHER" />
</intent-filter>
</activity>
<!-- 着信などを制御するBroadcastReceiver -->
<receiver android:name=".PhoneStateReceiver">
<intent-filter>
<action android:name="android.intent.action.PHONE_STATE" />
</intent-filter>
</receiver>

</application>

</manifest>

 ここで、前回のマニフェストに登場したBroadcastReceiverがやっと日の目を見ました。

 BroadcastReceiverを実装することで、以前ではアプリの起動時しか通話状態の変化に対応できなかったものが、バックグラウンド時でも対応できるようになったという訳でございます。

 次は、Toast表示を別の物に変更するということと、電話帳へのアクセスあたりをやりたいと思います。

PhoneStateListenerの作成 01

 android studio 2.3がリリースされました。

 Cameraアプリは難しいので、PhoneStateListenerを使ったアプリの作成に移ります。

 まだまだ未完成ですが、一歩進めたのでMainActivityを載せておきます。

MainActivity.java

import android.support.v7.app.AppCompatActivity;
import android.os.Bundle;
import android.telephony.PhoneStateListener;
import android.telephony.TelephonyManager;
import android.widget.Toast;

public class MainActivity extends AppCompatActivity {

@Override
protected void onCreate(Bundle savedInstanceState) {
super.onCreate(savedInstanceState);
setContentView(R.layout.activity_main);

// TelephonyManagerインスタンスを生成
TelephonyManager telephonyManager = (TelephonyManager)getSystemService(TELEPHONY_SERVICE);
telephonyManager.listen(phoneStateListener, PhoneStateListener.LISTEN_CALL_STATE);
}

public PhoneStateListener phoneStateListener = new PhoneStateListener() {
@Override
public void onCallStateChanged(int state, String incomingNumber) {
super.onCallStateChanged(state, incomingNumber);
phoneCallEvent(state, incomingNumber);
}
};

// 通話状態に応じて、Toastを表示する
private void phoneCallEvent(int state, String incomingNumber) {
switch (state) {
// 着信時の表示内容
case TelephonyManager.CALL_STATE_RINGING:
Toast.makeText(this, "着信中!" + incomingNumber, Toast.LENGTH_LONG).show();
break;
// 通話中の表示内容
case TelephonyManager.CALL_STATE_OFFHOOK:
Toast.makeText(this, "通話中!" + incomingNumber, Toast.LENGTH_LONG).show();
break;
}
}
}

 要するに、端末に着信が入った時と、端末が通話中の時にToast表示されるという単純な仕組みです。

 当然、通話状態などの端末情報にアクセスするために、マニフェストをこのように書き換えます。

AndroidManifest.xml

<?xml version="1.0" encoding="utf-8"?>
<manifest xmlns:android="http://schemas.android.com/apk/res/android"
package="com.example.ma2ri.telephonecontrol">

<!-- 電話番号、通話の状態などの端末情報を取得するためのパーミッション -->
<uses-permission android:name="android.permission.READ_PHONE_STATE" />

<application
android:allowBackup="true"
android:icon="@mipmap/ic_launcher"
android:label="@string/app_name"
android:supportsRtl="true"
android:theme="@style/AppTheme">
<activity android:name="com.example.ma2ri.telephonecontrol.MainActivity">
<intent-filter>
<action android:name="android.intent.action.MAIN" />

<category android:name="android.intent.category.LAUNCHER" />
</intent-filter>
</activity>
<!-- 着信などを制御するBroadcastReceiver
<receiver android:name="com.example.ma2ri.telephonecontrol.MainActivity$InComingCall">
<intent-filter>
<action android:name="android.intent.action.PHONE_STATE" />
</intent-filter>
</receiver> -->

</application>

</manifest>

 マニフェストの最後にBroadcastReceiverのコメントアウトがありますが、気にしないで下さい。まだまだ試行錯誤中です。

 ただ、この状態だと、待ち受けの状態にならないんです。そこを次に修正したいと思います。

時系列データ解析 03~ 復習

 少し駆け足でしたので、主に式の展開を詳しく説明します。

 まず、確率変数  x の期待値はこのように書けます。

 \displaystyle E ( x ) = \overline{ x } = \sum_{i = 1}^{n} x_{i} = \frac{x_{1} + x_{2} + \cdots + x_{n}}{n}

 分散はこのようになります。

 \displaystyle V ( x ) = E \left[ \left\{ x - E (x) \right\}^{2} \right] = \frac{(x_{1} - \overline{x}) + (x_{2} - \overline{x}) + \cdots + (x_{n} - \overline{x})}{n - 1}

 \displaystyle V ( x ) = E ( x^{2} ) - \left\{ E ( x ) \right\}^2

  i 番目の変数との偏差を  n - 1 あるいは  n で割ったものが分散です。分散の平方根標準偏差といいます。

 以上は  1 変数ですが、 2 変数の場合、分散は変数同士の関連具合を表す共分散となります。

 \displaystyle Cov ( x, y ) = E \left[ \left\{ x - E( x ) \right\} \left\{ y - E ( y ) \right\} \right] = E \left[ ( x - \overline{x} ) ( y - \overline{y} ) \right]

 \displaystyle Cov ( x, y ) = E( x y ) - E ( x ) E ( y )

 共分散は分散を一般化したもので、自己共分散は共分散の時系列データバージョンです。

 あらためて、 k 次の自己共分散  \gamma_{k} を記述します。

 \displaystyle \gamma_{k t} = Cov ( t_{t} , y_{t - k} ) = E \left[ ( y_{t} - \mu_{t} ) ( y_{t - k} - \mu_{t - k} ) \right]

 変数  y_{t} k 次の変数  y_{t - k}との関連具合を表しています。

 続いて、相関係数はこのように記述されます。

 \displaystyle \rho_{x y} = \frac{Cov (x, y)}{\sqrt{V(x)} \sqrt{V(y)}}

 相関係数(ピアソンの積率相関係数)は、共分散をそれぞれの標準偏差で基準化したものになります。もう少し詳しく書くとこうなります。

 \displaystyle \rho_{x y} = \frac{\sum ( x - \overline{x} ) ( y - \overline{y} )}{\sqrt{\sum( x - \overline{x} )^{2}}\sqrt{\sum ( y - \overline{y} )^{2}}}

 相関係数 \rho \leq 1 となりますが、それは上の式をコーシー・シュワルツの不等式( \mid (x, y ) \mid \leq \mid x \mid \cdot \mid y \mid)により展開することで確認できます。

 \displaystyle \mid \sum ( x - \overline{x} ) ( y - \overline{y} ) \mid \leq \sqrt{\sum( x - \overline{x} )^{2}}\sqrt{\sum ( y - \overline{y} )^{2}}

 \frac{\mid \sum ( x - \overline{x} ) ( y - \overline{y} ) \mid}{\sqrt{\sum( x - \overline{x} )^{2}}\sqrt{\sum ( y - \overline{y} )^{2}}} \leq 1

 ここで、 k 次の相関係数をあらためて記述します。

 \displaystyle \rho_{k t} = \frac{Cov ( y_{t}, y_{t - k} )} { \sqrt{V (y_{t}) \cdot V ( y_{t - k} )}}

 時系列データの自己相関係数相関係数と同じものであることが分かるのではないかと思います。