2018-04-01

隼時系列本：ノート5

本読み

以下の本を読みます。

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装
馬場真哉

プレアデス出版 2018-02-14
売り上げランキング : 7742

Amazonで詳しく見る by G-Tools

※ 以下、キャラクターが会話します。原作とは関係ありません。上の本からそれる話も多いです。誤りがあればご指摘ください。

前回：ノート4 ／次回：まだ

f:id:cookie-box:20180305231302p:plain:w60

モデルを選択したら残差の自己相関や正規性をチェックするって話だったけど、本に紹介されている他にはどんな方法があるんだろう？なんでこの方法なのかって気になるし…。

f:id:cookie-box:20180305232608p:plain:w60

ちょっと調べてみましょうか。まず自己相関の有無の検定から。主に Wikipedia を参考にしました。
上2つの検定は「 $1 \sim h$ 次の自己相関が全て0」を帰無仮説としていて、このような検定を総称してかばん検定（portmanteau test）というようですね。なぜかばんというのかはよくわからなかったんですが…。

Box–Pierce test （Ljung-Box test 原論文: https://www.researchgate.net/publication/246995234_On_a_Measure_of_Lack_of_Fit_in_Time_Series_Models）	帰無仮説が正しいならば、サンプルサイズ $n \to \infty$ の極限で $(\rho_1, \rho_2, \cdots , \rho_h)$ が平均ゼロ、分散 ${\rm Var}(\rho_k)=(n-k)/(n(n+2))$ 、共分散ゼロの多変量正規分布にしたがうことに着目すると、検定統計量を $Q=n \sum \hat{\rho}_k^2$ とすればこれが近似的にカイ2乗分布にしたがうので、この $Q$ が大きいときに帰無仮説を棄却しよう（何か自己相関があると疑おう）というものです。ただ、これは近似の仕方の関係で、自己相関があるにもかかわらず $Q$ が小さくなってしまう場合があったようなんです。これを改善したのが Ljung-Box test です。
Ljung–Box test （同上）	Box–Pierce test の検定統計量を $Q=n(n+2)\sum(n-k)^{-1} \hat{\rho}_k^2$ に修正したものです。本でメインに紹介されていますね。
Breusch–Godfrey test （Breusch–Godfrey test - Wikipedia）	これは時系列モデルの残差というより、線形回帰の残差に対して自己相関の有無を判定する検定のようなんですが…線形回帰の残差系列 $u_t$ を AR(p) でフィッティングしたとき、1～p次の自己相関が全てゼロならばこのフィッティングの決定係数がカイ2乗分布にしたがうことを利用するようです。
Durbin–Watson test （Durbin–Watson statistic - Wikipedia）	線形回帰の残差に1次の自己相関があるかどうかを検定します。この本では125ページで出てきますね。上の3つの検定が色々な次数の自己相関がゼロであることを一度に検定する一方で、1次の自己相関のみを検定というのは随分限定的なようですが、上の3つの検定が対応できない、誤差の分散が不均一な場合などにも対応できるようです。

次に正規性の検定です。これらの検定の帰無仮説は「データが正規分布にしたがう」ですね。正規性の検定は古くより数多く考案されており、初期には尖度や歪度が着目されていたようです。
Normality test - Wikipedia

D'Agostino's K-squared test	データの尖度と歪度がどれだけ正規分布のそれより逸脱しているかを $K^2$ 統計量としてまとめて検定するものです。 $K^2$ 統計量は近似的にカイ2乗分布にしたがいます。
Jarque–Bera test	本で紹介されていますね。考え方が D'Agostino's K-squared test と同じですが、検定統計量の近似がより洗練されているようです。
Anderson–Darling test （http://www.hep.caltech.edu/~fcp/statistics/hypothesisTest/PoissonConsistency/AndersonDarling1954.pdf）	この検定と Kolmogorov–Smirnov test は、正規性の検定というよりはもっと一般的に、データがある特定の確率分布から生成されたものかを検定します。もしデータがある確率分布 $f(x)$ から生成されているならば、そのデータのaパーセンタイル点と $f(x)$ のaパーセンタイル点は近いだろうというのはもっともなアイデアですよね。昇順にソート済のデータ $\{ x_i \}$ について、「そのデータがデータ中で小さい方から何%か」を横軸に、「その値は $f(x)$ の何パーセンタイル点か」を縦軸にプロットすると、 $(0, 0)$ から $(1, 1)$ まで直線上に並ぶだろうということです。この直線からのずれの2乗を足し合わせたものが Anderson–Darling test の検定統計量に相当しますが、帰無仮説の下でこのずれの分散は $F(x)(1-F(x))$ となり、分布の裾ほど小さくなります。なので、検定統計量は各点の誤差をこの分散の逆数で重み付けしています。この重み付けは、サンプリング誤差を均しているという見方もできますし、分布の裾のずれほど重くみているという見方もできます。 Anderson–Darling test によって正規性を検定するときは、平均／分散が既知／未知かでバリエーションがあることに注意ですね。
Kolmogorov–Smirnov test	Anderson–Darling test に先立つ手法ですね。この検定統計量は単に経験分布の累積分布関数と検定したい分布の累積分布関数の差のL∞-ノルムですね。なので、「2つのデータセットが同じ分布から生成されているか」という検定にもそのまま利用できます。2つの経験累積分布の差に置き換えるだけですので。これは画一的で便利な検定である一方、正規性の検定としての検定力は Anderson–Darling test や Shapiro-Wilk test に劣るそうですので、これらの検定が利用できる場面ではそちらを利用した方がよさそうですね。
Shapiro-Wilk test	これは分布全体の形というよりもそれぞれの値の位置に着目していますね。サンプルサイズが $n$ だったとき、小さい方から $i$ 番目のデータが、正規分布から生成された $n$ 個のサンプルの $i$ 番目に小さいデータの期待値に近いか、といったイメージです。2011年のある研究によると、正規性の検定の中で最もよい検出力をもつとか。検出力がよいということはつまり、データが正規分布から生成されたものでないときに、ちゃんと帰無仮説を棄却できる確率が高いということですね。一方で、データの中にタイ（複数のデータが全く同じ値となる）が多い場合に弱いなど、弱点もあります。小さい方から $i$ 番目のデータの位置、などに着目しているので、タイには弱そうですよね。

それを聞くと自己相関の有無の検定はやっぱり Ljung–Box test がスタンダードなのか。正規性の検定もたくさんあるな。同じことを確認しようとするのに、色々な考え方があるってのが面白いけど面倒かも。

統計的検定では、帰無仮説を棄却できるかどうか、つまり、目の前のデータがあるパラメータをもつ確率分布から生成されていることを否定できるかどうかを判定しようとします。「帰無仮説が正しくなさそうな度合い」を何で測るかによって色々な考え方が生じるでしょう。

でも目の前のデータの背後に確率分布があるっていうのが人間の空想なのに、この確率分布じゃないと否定されても、確率分布にしてみればとばっちりだよな。

え？すみませんが、棄却された確率分布の側がどう感じられているかは僕にはちょっと…。

72ページからは時系列分析を実際にやってみようって話か…といっても今パソコンないしな…。

ノートPCならありますよ。

なんで学校にノートPC持ってきてんの！？

ちょっとした編曲作業に使おうと思いまして。学校でアイデアが浮かぶこともあるかもしれませんし。

いや確かに曲のアイデアは逃がしたくないけどさ…。

でもこれ再現実行しても本の内容なぞるだけになっちゃいますし、とばして第3部を読みましょうか。見せかけの回帰を別の具体例でやってみましょう。

そうそう、第3部に「見せかけの回帰」ってあったけど、何が見せかけなんだ？

見てみるのが早いかと。互いに全く関係ない2つの時系列モデル

$x_t = 0.8 x_{t-1} + 0.1 x_{t-2} + 0.2 \varepsilon _{x, t-1} + \varepsilon _{x, t}$
$y_t = 0.6 y_{t-1} + 0.2 y_{t-2} + 0.1 \varepsilon _{y, t-1} + \varepsilon _{y, t}$

を用意しましょう。たまたまどちらも ARMA(2,1) にしてみました。では、これらのモデルから 100 点ずつデータを生成し $(x_1, \, y_1) , \cdots, (x_{100}, \, y_{100})$ を得て、線形回帰するとどうなるでしょう。

どうなるって…互いに全く関係ないデータに線形回帰しても意味ないだろ。

実際にやってみましょう。データを生成して R の lm() で線形回帰分析します。lm() はデータに線型モデルを仮定すると、モデル係数の最小2乗推定量や残差、決定係数やF値などを計算してくれる関数ですね。以下では $y_t = \alpha + \beta x_t + \varepsilon_t \; \bigl( \, \varepsilon_t$ は互いに独立に $N(0, \sigma^2)$ にしたがう $\bigr)$ を仮定することになります。

set.seed(1234)
x <- arima.sim(n=100, model=list(ar=c(0.8, 0.1), ma=c(0.2)))
y <- arima.sim(n=100, model=list(ar=c(0.6, 0.2), ma=c(0.1)))
plot(as.numeric(x), as.numeric(y))
print(summary(lm(y~x)))

実行結果は以下です。傾きが約 -0.4、切片が約 1 となりました。傾きのp値も切片のp値も小さく、傾きも切片も有意であるという結果になっています。つまり、傾きも切片もゼロとは思えないということですね。

Call:
lm(formula = y ~ x)

Residuals:
     Min      1Q  Median      3Q     Max 
 -3.7131 -0.9341 -0.1749  0.9870  3.7385 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.06548    0.16071   6.630 1.84e-09 ***
x           -0.43461    0.08694  -4.999 2.52e-06 ***

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.454 on 98 degrees of freedom
Multiple R-squared:  0.2032,    Adjusted R-squared:  0.1951 
F-statistic: 24.99 on 1 and 98 DF,  p-value: 2.523e-06

え…？ $x$ も $y$ も互いに関係ないモデルから生成されたのに、線型モデルをあてはめると傾きや切片がゼロとは思えないって…どうして？？

これが見せかけの回帰ですね。

えー何が起こってるのか全然わかんない！これって、p値的にはそのモデルで表される線形関係があるはずって結果なんだよな？じゃあもう見せかけの回帰じゃなくてリアル回帰じゃん。

まだ気付きませんか、ハヤト。バンドのリーダーなんだからしっかりしてください。

バンド関係ない！

線形回帰分析の係数のt検定の検定統計量がt分布にしたがうのは、あくまで「残差が互いに独立で平均ゼロの正規分布にしたがう」という仮定の下での話です。 $x_t$ と $y_t$ が自己回帰モデルから生成される場合にこの仮定は成り立つでしょうか。簡単のために、どちらもランダムウォークとして計算してみます。

$x_t = x_{t-1} + \varepsilon_{x, t}$
$y_t = y_{t-1} + \varepsilon_{y, t}$
$y_t = \alpha + \beta x_t + \varepsilon_t \; \Rightarrow \; \varepsilon_t = y_t - \alpha - \beta x_t$
$\varepsilon_t = y_{t-1} + \varepsilon_{y, t} - \alpha - \beta (x_{t-1} + \varepsilon_{x, t}) = \varepsilon_{t-1} + \varepsilon_{y, t} - \beta \varepsilon_{x, t}$

これより、残差 $\varepsilon_t$ もランダムウォークになることがわかります。残差がランダムウォークするのに、最小2乗推定で回帰係数を推定するべくもありません。ましてや、検定統計量がt分布にしたがうはずがないんです。これを考慮せずに lm() で線形回帰分析しても、その結果はまやかしであり、紛れもなく「見せかけ」の回帰なんです。

t検定が利用できる状況が整っていないのにt検定をしたのがいけなかったのか。じゃあどうすればいいんだ？

残差に自己相関があるのにt検定をしてしまうのが悪いんですから、残差に自己相関が残らないように階差系列をとるなどデータの方を加工するか、回帰モデルの方で自己相関を表現できるように修正するか、あるいは残差に自己相関がある可能性を踏まえた回帰分析の方法を導入する必要があるでしょう。先ほどのデータに一般化最小2乗法（GLS）を適用してみましょう。まず本で紹介されている prais パッケージの関数から。

library(prais)
print(prais.winsten(y~x, data.frame(x=x, y=y), iter=1))

結果はこうです。もはや傾きのp値が小さくありませんね。

Call:
lm(formula = fo)

Residuals:
     Min      1Q  Median      3Q     Max 
 -3.3421 -0.5651  0.0661  0.5300  3.0788 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)  
Intercept  0.76130    0.32288   2.358   0.0204 *
x         -0.04359    0.10781  -0.404   0.6869  

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9879 on 98 degrees of freedom
Multiple R-squared:  0.05423,   Adjusted R-squared:  0.03493 
F-statistic:  2.81 on 2 and 98 DF,  p-value: 0.06508

      Rho Rho.t.statistic Iterations
 0.688919         9.44636          1

他にも GLS は色々なパッケージがあるようなので、少し調べて出てきたのもやってみます。rms パッケージの Gls() ですが、どのような自己相関をもつのか引数 correlation に明示的に与えるようです。この引数を省略すると回帰は有意となってしまいますが、手抜きで corAR1() と指定してみると、回帰は有意ではなくなりました。 Gls function | R Documentation

library(rms)
print(Gls(y~x, data.frame(x=x, y=y), correlation=corAR1()))

次のようなページも見つけました。
一般化最小二乗法と識別問題 - hnami.net_Pukiwiki

ふーん、とりあえず残差に自己相関がある場合の回帰にも対処できる方法があるならよかった。…138ページからの、共和分って何？

簡単な例だと、 $x_t$ と $y_t$ がそれぞれ単位根過程だったとして、 $x_t$ と $y_t$ を適当に線形結合した過程が定常になるなら、 $x_t$ と $y_t$ は共和分しているといいます。

なるほど…だから何？

単位根過程は差分系列をとって定常過程にしてから解析するのがセオリーです。しかし、 $x_t$ と $y_t$ が共和分している場合、先にそれぞれ差分系列をとってしまうと、 $x_t$ と $y_t$ の関係が失われてしまうということですね。下の沖本氏の本の132ページから、共和分を踏まえなければ導出できない経済事象の具体例が3つ挙げられています。
経済・ファイナンスデータの計量時系列分析 (統計ライブラリー) | 沖本竜義 |本 | 通販 | Amazon
つまり、1つ目の例は消費系列と所得系列の関係を調べるのに、これらが共和分していることそのものが示唆になっているのですが、先に差分系列をとってしまうとこの関係が見えなくなってしまいます。共和分している系列どうしに回帰モデルをあてはめると、残差が定常になりますので、共和分が疑われるときはそれを検定することになります。本の142ページをみると、この用途専用の検定をつかわなければならないようですね。

（ノート6があれば）つづく

2018-03-13

隼時系列本：ノート4

本読み

以下の本を読みます。

※ 以下、キャラクターが会話します。原作とは関係ありません。上の本からそれる話も多いです。誤りがあればご指摘ください。

前回：ノート3 ／次回：まだ

時系列のモデルを決めるのに、次数は AIC か何かでゴリゴリ決めるけど、何階差分系列を取るべきかはそれじゃ駄目で、仮説検定をつかって対象データが単位根かどうか判断するって話で…KPSS検定とADF検定は帰無仮説と対立仮説が逆？

有名な単位根検定の中でKPSS検定だけが「単位根がない」の側を帰無仮説とするようですね。なるほど「単位根がある」の方を帰無仮説とするのが自然なように思います。両側検定では $H_0 : \mu = \mu_0 \, \; {\rm vs.} \, \; H_1: \mu \neq \mu_0$ の形式で $H_0$ が正しいかどうかを確かめますが、「単位根がある」とはつまり「その時系列を生成する確率過程を AR(1) で表現したとき $\phi_1 = 1$ である」ですので、 $H_0: \phi_1 = 1$ になりそうです。しかし、「単位根がある」を帰無仮説とするDF検定やADF検定では、帰無仮説を棄却できないケースに困ったようなんです。単位根がないというよほど「強い証拠」がない限り帰無仮説が受容されてしまうと（以下；KPSS検定の原論文のイントロダクションより）。

Testing the null hypothesis of stationarity against the alternative of a unit root

「単位根があることを棄却できない」ってことは「単位根がないと結論づけるには至らない」ってことだよな…それは困るな。棄却されたら単位根過程じゃないって安心できるけど、棄却されなかったらどっちなんだよっていう。

なので、KPSS検定の提唱者らは、帰無仮説と対立仮説を逆にしたんです。彼らは最初からランダムォーク項を仮定し、「ランダムウォークの分散がゼロである」を帰無仮説としました…本の64ページにある通りです。

なるほど、それなら帰無仮説が棄却されれば「単位根がある」ことになる…別の検定では単位根があるのかないのかどっちかわからなかったケースでも、きっぱりと「単位根がある」っていえるようになるかも。…じゃあKPSS検定だけでよくない？何階差分をとるか決めるときには「単位根がある」ってはっきり言ってくれる検定の方がうれしいし…。

実際この本ではKPSS検定を時系列モデル構築手順に据えていますが、65ページにデータに合わせて検定方法を変えるようにとありますね。どの検定を採用すべきかはどのような仮定を置いているかに依存しますし、時系列データの長さや、 $\phi_1$ の大きさにも依存するとも（ $\phi_1 = 0.1$ ならなかなか単位根があるようにはみえないかもしれませんが $\phi_1 = 0.9$ なら単位根があるようにみえてしまいそうですよね）。以下のような論文をWeb上で見つけました。

https://www.czso.cz/documents/10180/32912822/32019716q3047.pdf/09710b90-e1d0-4bb1-816e-5b83faad686b?version=1.0

結局どの検定にするかもケースバイケースってことか…そういえば、65ページのDF検定に「簡単そうに見えますが（略）うまくいきません」ってあったけど、まず簡単そうに見えるって気持ちがわかんない！

DF検定って式の形だけみると単なる単回帰分析にみえます。単回帰分析とは… $(x, y)$ のデータ対がたくさんあったとします。(気温, アイスクリームの売上高) のデータと考えてもいいです。このデータから傾き「1度気温が上がると、アイスの売上高がどれだけ伸びるか」を求めたいとします。単純な仮定の置き方の1つは、 $y$ が確率分布 $N(ax+b, \, \sigma^2)$ から生成されると決めることです。最小2乗法で $a, \, b$ の推定値 $a', b'$ を求めることができるでしょう。全てのデータが $xy$ 平面で $y$ 軸に平行に並んでいるのでもない限り、 $a'$ の値は必ず何かしら求まります。しかし、 $a'$ が求まるからといって先の仮定が正しいとは限りませんよね。本当は直線に当てはめるべきデータではないかもしれません。なので、 $H_0: a=0$ の下でその $a'$ が得られるかどうか確かめるんです。もし $H_0$ が棄却されれば、傾きが $a'$ であることが信頼できるでしょうし、棄却されなければ最初の仮定が信頼しづらいということです。…とまあここまでは簡単な話ですね。ただ、DF検定は単回帰分析のようにみえてこのように簡単にはいかないということなんです。

ごめんその話も別に簡単じゃないんだけど…何がさらに難しいの…。

$y_t= a x_t + \varepsilon_t$ と $y_{t} = \phi_1 y_{t-1} + \varepsilon_t$ は違うんです。後者では、 $\varepsilon_{t=0}$ の影響がいつまで経ってもなくなりません。ので、 $a$ の推定値 $a'$ の分布と、 $\phi_1$ の推定値 $\phi_1'$ の分布の形は違ってきます。さらに、いま帰無仮説は $H_0: \phi_1 = 1$ ですから、この非定常な場合での分布の形を知らないといけないんです。

…よくわかんないけど、単位根過程かどうか調べるのに苦労してるんだな。しかもそれだけ苦労したのにADF検定は帰無仮説が棄却できなかったりするんだろ？

なので、改良されたADF-GLS検定というのもあるそうです。と以下の資料に。

https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/155556/1/kronso_183_2_35.pdf

へー。まあそれでARIMAモデルが決まったとして、残差をチェックするのか。…残差に自己相関が残ってたら駄目っていうけどさ、計算機が決めてくれたのに ARMA モデルの次数が最適になっていないってあるの？それとも残差が残っちゃったらもう ARMA モデルの限界ってこと？

僕もよくわかりませんが…ARMAのパラメータ最適化も不安定な場合もあるのかもしれませんが…あるいは、階差を取るべきなのに取っていなかったなどということもあるかもしれませんね。

（ノート5があれば）つづく

2018-03-08

隼時系列本：ノート3

本読み

以下の本を読みます。

※ 以下、キャラクターが会話します。原作とは関係ありません。上の本からそれる話も多いです。誤りがあればご指摘ください。

前回：ノート2 ／次回：まだ

第2部の4章までで色々なモデルを扱って、第2部の5章からは、じゃあどのモデルを選べばいいかって話か…。

ここでモデルという言葉は ARIMA(2,0,0) と ARIMA(1,1,0) は異なるモデルということですね。なので、モデル選択というのは次数の決定のことであると書いてありますね。

あれ、AR モデルか MA モデルか ARIMA モデルか SARIMA モデルか、っていう選択を最初にするのかと思ったんだけど、その選択はないの？

まず、ARIMA モデルは AR モデルと MA モデルを包含しますよ。SARIMA モデルはさらに ARIMA モデルを包含しますが、周期を与えてつかうモデルなので、対象データに周期がある場合には SARIMA をつかってみればいいんじゃないですか。

あ、そうか。じゃあやっぱり最初に取りかかるのは次数の選択なのか。確かにARIMAモデルをつかうとしても次数をどうすればいいのか全然わからないよな…これはどういう手順で決めるのか気になるぞ…えっ、「手当たり次第」？ショックなんだけど…。

ショックを受けることないでしょう。僕たちは計算機が発達した時代に生まれたんです。考えなくても解決できる手順を中途半端に考える必要はありません。下手の考え休むに似たりです。しかし、計算機にまかせるにせよ、モデルのよさの指標については把握しておく必要があるでしょう。

う、うん…まずAICっていう指標が紹介されてるよ…それでその前に、尤度？尤度って何？この61ページのコイン投げの 2/9 っていう尤度は大きいの？小さいの？

コインの表が出る確率が 1/3 ではなく 1/2 だったらその尤度はどうなります？

それは、1/2×1/2＝1/4 になる。さっきの 2/9 と比べると 1/4 の方が大きい？

なので、あるコインについて「2回投げたら1回目が表で2回目が裏が出た」ということがわかっているとき、このコインの表が出る確率は 1/3 よりも 1/2 の方が尤もらしいということです。もっと一般に、コインの表が出る確率を $p$ とすると、「2回投げたら1回目が表で2回目が裏が出た」という結果に対する $p$ の尤度は $f(p)＝p(1-p)$ ですから、これを最大にするのは $p＝1/2$ です。なので結局、2回投げた結果しかわからない範囲では、このコインの表が出る確率として最も尤もらしいのは 1/2 です。

…2回投げて1回目が表で2回目が裏だったんだから、コインの表が出る確率が 1/2 らしいって当たり前なんじゃ？

まあ平均値については、「標本平均の期待値」が「真の平均値」に一致しますからね。でも平均値ではない分布のパラメータについては決して当たり前では…この話は脇道にそれるので置いておきましょう。…ところでハヤト、このコイン投げの結果から考えられる、もっと尤もらしいコインはあるでしょうか。

もっと尤もらしい？「2回投げたら1回目が表で2回目が裏が出た」という結果から考えられる？さっき尤度 $f(p)$ を最大にするのは $p=1/2$ だってジュン自分で言ってたじゃん。だから「表が出る確率が 1/2 のコイン」より尤もらしいコインはないだろ？

$f(p)$ はコインの表が出る確率が1回目も2回目も同じ場合の尤度ですので。2回目に投げるときは1回目に投げるときとはコインの性質が変化しているかもしれません。コインを特徴付けるパラメータを $p_1$ と $p_2$ の2個に増やして、1回目に投げるときは確率1で表が出て、2回目に投げるときは確率1で裏が出るコインを考えるとどうでしょう。こうすれば尤度は 1×1＝1 です。テキストの61ページと同様の例ですが。

えーそれはズルイだろ。

ええ、ずるいんです。ずるいというよりは、このモデリングは3回目に投げた結果を予測するのに全く有用ではないでしょう。なので、なるべくパラメータを増やさずに尤度を最大化したいんです。なるべくパラメータを増やさずに尤度を大きくできたかどうかの指標が62ページの AIC です。

62ページ…これがAICの式か…最大化対数尤度？

61ページに説明がありましたが、あるモデル、例えば ARMA(1,1) でも何でもいいです、ARMA(1,1) だと $c, \, \phi_1, \, \theta_1$ の3つのパラメータが推定対象ですね。対象時系列データの下での尤度が最大になるようにこれらのパラメータを最適化した上で、その尤度に対数をとったものが最大化対数尤度です。最後に対数をとる必要もありませんので、対数尤度が最大になるように3つのパラメータを最適化したときに対数尤度、と言った方がいいでしょうか。次数の異なるモデル間でよさを比較するのだから、パラメータセットはベストにしておかなければならないでしょう。

…時系列データに対するARIMAモデルの尤度ってどう計算するの？時系列データって、時点ごとにデータたくさんあるよな？

さっきのコイン投げの例でも、あるパラメータの下での2回のコイン投げの結果の同時確率がそのパラメータの尤度だったでしょう。2回だろうと100回だろうと一緒ですよ。与えられた100日分の気温に対するあるパラメータセットの尤度は $p($ 1日目～100日目の気温 $|$ パラメータセット $)$ です。もっとも、ARIMAモデルの尤度関数は…インターネットを探せばありますよ。例えばこのリンクなどがそうですね。

http://econ.nsysu.edu.tw/ezfiles/124/1124/img/Chapter17_MaximumLikelihoodEstimation.pdf

…まあいいや。その最大化対数尤度と、パラメータ数で AIC が決まるのか。尤度は大きいほどよくて、パラメータ数は小さいほどよくて、AICは小さいほどよい指標なんだな。…ジュンさっき、なるべくパラメータを増やさずに尤度を大きくしたいって言ってたよな。(2-41) 式をみるとさ、パラメータを1つ増やすんだったら、対数尤度は1よりも大きく増えてくれなきゃ困るってことだよな。なんでパラメータ1個と対数尤度1単位が釣り合うんだろう。あとなんで第1項も第2項も2がかかってるの？

以下のリンクに AIC の導出がありますよ。もっとも、AIC は原論文も探せばありますが。自由度 k のカイ2乗分布の期待値が k なので AIC は (2-41) 式のようになっているんですね。全体に2がかかっているのは、こうすると第1項が漸近的にカイ2乗分布にしたがうからですね。

http://www.is.titech.ac.jp/~shimo/pub/koza2000/ohp-20000703b.pdf

ふーん…まあなんやかんやで計算機が ARMA の次数を決めてくれるんだな。でも階差についてはそういうわけにいかないのか…それが63ページからの検定…帰無仮説って何？

仮説検定ですか。例えば有意水準5%で仮説検定するのだったら、「帰無仮説が正しいとしたら起きる確率が5%以下になるくらい珍しいこと」を「偶然それが起きたと考えるには低い」とみなすんです。与えられたデータが、「帰無仮説が正しいとしたらこんなデータが得られる確率は5%以下だ」というデータだったら、帰無仮説を棄却します。つまり、帰無仮説は正しくないと考えるんです。もし帰無仮説が正しいときに与えられたデータが得られる確率が5%よりも大きいのだったら、帰無仮説を受容します。これは、帰無仮説は正しくないとまではいえない、という状態ですね。

うーん…なんか具体例ない？

今調べたんですけど、高校2年生男子の平均体重は60.4kg、標準偏差は9.93kgらしいです。この平均体重が未知だったとして、僕とハヤトが高校2年生男子の無作為標本だったとき、帰無仮説「高校2年生男子の平均体重は60kg以上」は棄却されるでしょうか。有意水準5%とします。…分散既知の場合の平均の検定では、標本平均が正規分布にしたがうことを利用します。今回は標本平均が下側5%点より小さい領域が棄却域です。正規分布の下側5%点は -1.64 で、サンプル数2なので、棄却域は $60 - 1.64 \times \sqrt{9.93^2 / 2} = 48.48$ 未満ですね。僕のハヤトの体重の平均は 48.5kg なので棄却域になく、帰無仮説は棄却されません。帰無仮説が正しい下で5%も起こりえない珍しいデータではないということですね。

でもすごいギリギリ…。

（ノート4があれば）つづく

2018-03-07

隼時系列本：ノート2

本読み

以下の本を読みます。

※ 以下、キャラクターが会話します。原作とは関係ありません。上の本からそれる話も多いです。誤りがあればご指摘ください。

前回：ノート1 ／次回：まだ

えーと、前回やったのは、時系列分析っていうのは時系列モデルを推定することで、時系列データには周期性とかの構造があって、第2部が具体的なモデル推定の手順か（前回の最後にジュンがちょっとしゃべっちゃってたけど…）。最初はまず定常過程の説明からで、平均と分散が時刻によらず一定で、自己相関も時点によらず時間差にのみ依存するのが定常過程で、そうじゃないのが非定常過程ね。でも非定常過程でも定常過程に変換できることがあるんだ。変換できる例が単位根過程？なんか変な名前…。

単位根過程は1階差分が定常過程になるような非定常過程のことですね。そのような名前である理由まではこの本には触れられていないようですが、自己相関を特徴付けるある方程式が1を解にもつことに由来していて…その方程式については後の方で出てくるようですね。

ふーん？ …なあジュン、38ページの (2-6) 式って「定数＋ホワイトノイズ」なの？「定数＋ホワイトノイズ－1時点前のホワイトノイズ」にみえるんだけど、「ホワイトノイズ－1時点前のホワイトノイズ」は「ホワイトノイズ」ってこと？

確認すればいいでしょう。「ホワイトノイズ－1時点前のホワイトノイズ」の平均が0なのは自明、分散は ${\rm Var}(\varepsilon_t - \varepsilon_{t-1}) = {\rm E} \bigl( (\varepsilon_t - \varepsilon_{t-1})^2\bigr) = {\rm E}(\varepsilon_{t}^2) + {\rm E}(\varepsilon_{t-1} ^2)= 2 \sigma^2$ なので時刻によらず一定ですね。1次の自己共分散は ${\rm Cov}(\varepsilon_t - \varepsilon_{t-1}, \varepsilon_{t-1} - \varepsilon_{t-2}) = {\rm E} \bigl( (\varepsilon_t - \varepsilon_{t-1})(\varepsilon_{t-1} -\varepsilon_{t-2})\bigr) = {\rm E}(\varepsilon_{t-1}^2 )=\sigma^2$ なので1次のみゼロではありませんが、時点によらず時間差のみに依存するので結局定常過程ですね。ならいいじゃないですか。

でもちょっと気になるじゃん。サポートページにコメントしといた…ら返信もらえた！

何やってるんですか…。

42ページからはARIMAモデル？まずはARモデル…AR(1) だったらある時点の値は1時点前の値で回帰されるってことか。これって、昨日の気温が暑かったから今日も暑い、みたいな表現ができるってことだよな。iid 系列よりこういうモデルの方が時系列モデルって感じがするかも。…あれ？ iid 系列は定常過程だったけど、ARモデルは定常過程なの？ Box-Jenkins 法ではまず定常過程にしてから色々なモデルを適用するんだよな。このARモデルもそんな風にしてから適用するモデルの例ってことは、定常過程ってこと？

44ページに、AR(1) の特別な場合について何て書いてあります？

えっと待って…係数が1のときランダムウォーク？ランダムウォークって AR(1) の特別な場合なんだ。ランダムウォークは非定常過程だから、あれ、ARモデルは非定常過程ってこと？？

とは限りません。一般の AR(1) の期待値、分散、自己共分散を考えてみましょう。AR(1) の一般形は $y_t = c + \phi_1 y_{t-1} + \varepsilon_t$ です。両辺の期待値をとると、 ${\rm E}(y_t) = c + \phi_1 {\rm E}(y_{t-1})$ になり、両辺の分散を取ると ${\rm Var}(y_t) = \phi_1 ^2 {\rm Var}(y_{t-1}) + \sigma ^2$ になりますね。もしこれが定常過程であれば、 ${\rm E}(y_t) = {\rm E}(y_{t-1}) = \mu, \; {\rm Var}(y_t) = {\rm Var}(y_{t-1}) = \gamma$ ですが、 $|\phi_1| < 1$ のとき $\mu = c / (1 - \phi_1), \; \gamma =\sigma^2 / (1 - \phi_1^2)$ はこれを満たしますね。このとき自己共分散はどうでしょうか。 $k$ 次自己共分散は ${\rm Cor}(y_t, y_{t-k}) = {\rm Cor}(\phi_1 y_{t-1} + \varepsilon_t, y_{t-k}) = \phi_1 {\rm Cor}(y_{t-1}, y_{t-k}) = \phi_1^2 {\rm Cor}(y_{t-2}, y_{t-k}) = \cdots = \phi_1 ^k \gamma$ となります。これは時点に依存しませんので、定常過程です。つまり、AR(1) は $|\phi_1| < 1$ ならば定常過程になりえます。逆に $|\phi_1| < 1$ ならば期待値や分散は上の値に収束し、結局必ず定常過程になるのですが。

ARモデルが定常か非定常かは係数によるってことか…なんか面倒だな。定常データのモデリングにつかうつもりなら、もう $|\phi_1| < 1$ ってことに最初から決めておけばよくない？

1次のARモデル AR(1) の場合はそれでよいんですが…2次以上になると係数への制約をそう簡単に書き下せないんです。先ほどの AR(1) の $k$ 次自己共分散は、 $k-1$ 次自己共分散と ${\rm Cor}(y_t, y_{t-k}) = \phi_1 {\rm Cor}(y_{t-1}, y_{t-k})$ という関係が成り立っていましたよね。もっと一般に、p次のARモデル AR(p) の $k$ 次自己相関 $\rho_k$ について、 $\rho_k = \phi_1 \rho_{k-1} + \phi_2 \rho_{k-2} + \cdots + \phi_p \rho_{k-p}$ が成り立つんです。この式は、 $\rho_k$ からはじめて順に $\rho_{k-1}, \, \rho_{k-2}, \, \cdots, \rho_{0}=1$ まで自己相関を紡いでいく漸化式と見做すことができるでしょう。この漸化式の特性多項式は $f(z) = 1- \phi_1 z - \cdots \phi_p z^p$ です。同じ式が本の50ページにありますね。 $\rho_k$ はこの特性多項式の $p$ 個の根を $-k$ 乗したものの線形和で書けるんです。漸化式は数学Bでやりましたよね？数列 $a_n$ の隣接3項間漸化式を特性多項式の根を用いて解くときに、 $a_n$ の一般形に根 $\alpha$ の $n$ 乗があらわれたのを思い出してください。ここでは（添え字の向きがややこしいですが） $\rho_{0} \propto \alpha^k \rho_{k}$ かつ $\rho_0 = 1$ なので $\rho_k \propto \alpha^{-k}$ ということです。ということは、もし根 $\alpha$ の中に絶対値が1以下のものが混じっていたら、 $k$ が大きくなるほど $\alpha ^{-k}$ が振動か爆発し、ARモデルは不安定になってしまいます。というわけで、特性多項式のすべての根の絶対値が1より大きいことが、ARモデルが定常である条件なんです。

…よくわかんないけど、ARモデルが定常になる条件が面倒なのはわかったよ…。

ちなみに、37ページで出てきた「単位根過程」の名前はこのARモデルの特性多項式の「根」に由来します。単位根過程はARモデルで表現したとき特性多項式の根の1つが1になるんです。これは定常過程の1階和分を取って特性多項式を確かめてみるといいと思います。例えば―

あー！それは置いとこう！なんか全然話進んでないし！本に戻ると、ARモデルとは別にMAモデルっていうのもあるって。これは…過去 q 時点のノイズで現時点の値を回帰するってことか。…でも、ARモデルも過去の値を通して過去のノイズが含まれているよな。どう違うんだろう？

…すみません、確かに脇道にそれ過ぎましたね。ARモデルとMAモデルの関係は49ページにありますね。AR(1) は実質 MA(∞) だと。だから、過去の実現値そのものと関係がありそうな時系列データなら、ARモデルが向いているでしょう。無理やり MA(q) モデルで推定すると、q を非常に大きくしなければならないでしょうから。逆にMAモデルをARモデルで表現するときは…この話は50ページにありますね。MAモデルは（同じ期待値、分散、自己相関をもつ過程が複数存在するのですが、その内「反転可能」といわれる1つの過程は）AR(∞) に書き直すことができます。なので、MAモデルから生成されたデータを無理やり AR(p) モデルで推定しても p が大きくなってしまい、よくないでしょう。ARモデルとMAモデルはどちらかがどちらかを兼ねるというわけにはなっていないんですね。

なるほど。47ページにはARモデルとMAモデルのコレログラムがある…これは、これ見てどう思えばいいの？

モデルを推定するのに、対象時系列データの自己相関や偏自己相関のコレログラムを眺めてモデル選択するのが伝統的な手順だったんです。しかしこの本は、そう教えていませんね。そのようなコレログラムの観察は、計算資源が乏しかった頃の「古い同定の手順」だと60ページに。ただ、伝統的な手順にならえば…例えばハヤトが分析したい対象時系列データの標本自己相関をコレログラムにプロットしたところ、2次以降の自己相関がほとんどゼロだったとします。ARモデルとMAモデルのどちらを適用したいですか？

自己相関？って偏自己相関じゃない方だよな？ 47ページだと ACF っていうグラフの方…がすぐにゼロになっているのは…MAモデルの方？

ええ、逆に偏自己相関の方がすぐにゼロになっていたらARモデルを適用、という具合です。MA(q) は q+1 次以降自己相関をもちません。AR(p) は p+1 次以降偏自己相関をもちません。しかし、MA(q) は AR(∞) なので無限に偏自己相関をもち、AR(p) は MA(∞) なので無限の自己相関をもつんですね。

ARモデルとMAモデルってそういう違いがあるのか。それで、両方を組み合わせたのがARMAモデルで、ARIMAモデルというのは、何度か差分を取ってからARMAモデルを適用するってことか。54ページからのSARIMAモデルは、季節成分を考慮したARMAモデル…え、(2-39) 式って何これどうなってんの？？

落ち着いてください。例えば SARIMA(0,0,1)(0,0,1)[12] を考えてみましょう。このとき、(2-39) 式は、 $y_t = \theta(B) \Theta(B) \varepsilon_t = (1+\theta_1 B)(1+\Theta_1 B^{12}) \varepsilon_t = (1 + \theta_1 B + \Theta_1 B^{12} + \theta_1 \Theta_1 B^{13} ) \varepsilon_t$ ですから、ラグ演算子 $B$ を用いずに書けば、 $y_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \Theta_1 \varepsilon_{t-12} + \theta_1 \Theta_1 \varepsilon_{t-13}$ ということです。

なんだそういうことか…じゃあ、SARIMA(0,0,1)(0,0,1)[12] は、1時点前のノイズに依存していてその係数が $\theta_1$ で、12時点前のノイズにも依存していてその係数が $\Theta_1$ で、13時点前のノイズにも依存していてその係数が $\theta_1 \Theta_1$ なのか…。なんか13時点前って気持ち悪くない？「1時点前」はさっきのノイズへの依存って感じで、「12時点前」は昨年のこの月のノイズへの依存って感じでわかりやすいけど、「13時点前」ってさあ…。

昨年のこの月の影響を織り込み済みのノイズ「 $\varepsilon_t + \Theta_1 \varepsilon_{t-12}$ 」を新たに「 $\tilde{\varepsilon_t}$ 」とみなしてこれで MA(1) を適用するって考えればいいんじゃないでしょうか。

あー確かに。SARIMAモデルって（月ごとに値がある周期12の時系列だったら）こんな感じなのかな？

必要ならデータ全体を前年同期との（D階）差分にする。
必要ならデータ全体を前月との（d階）差分にする。
同じ月の系列「2000年1月の値、2001年1月の値、2002年1月の値、…」を説明する ARMA(P, Q) モデルを決める（但し、どの月も同一のモデル）。
さっきの ARMA(P, Q) での過去のノイズへの依存も含めたノイズを新しい「現時点のノイズ」とみなし、過去の値への依存も含めた値を新しい「現時点の値」とみなし、ARMA(p, q) を適用する。

でもこれってなんかパラメータ多いけど、同じ期待値や自己相関をもつ過程が1つに決まるのかな…。

SARIMA が実装されているパッケージでは上手くやっているんじゃないでしょうか…。

第2部の4章の最後は、ARIMAXモデル…これはイベントの影響を考慮することができるモデルで…ダミー変数って何？過去の値とか過去のノイズだけじゃなくて、こういうのにも依存していいの？

そうですね、 $y_t$ が各日のドーナツ屋の売上高だったとしたら、よりよく推定するために、売上高以外にも日ごとの情報を色々追加してしまおうということですね。 $x_{k, t}$ は「その日が平日か休日か」「その日が晴れか雨か」「その日に近くでコンサートがあるか」などがありえると思います。しかし、「平日」も「休日」も数字ではありません。そこで仮に「平日＝0」「休日＝1」などと適当な数字を割り当ててしまいます。これがダミー変数です。

（ノート3があれば）つづく

2018-03-05

隼時系列本：ノート1

本読み

以下の本を読みます。

※ 以下、キャラクターが会話します。原作とは関係ありません。内容の誤りは本ブログ筆者に帰属します。

ジュン、忙しいところ悪いけど、この本一緒に読んでくれない？

何ですかハヤト…時系列分析？何だってそんな本読む気になったんです？

これ先月出た本なんだけど、表紙の絵がハヤブサで、通称「隼時系列本」っていうらしいんだ。それでプロデューサーが、「隼人が隼時系列本を読む」っていうネタを思いついちゃって…。

導入が雑すぎる…。

…それで読もうとしたんだけど、まずこれは何の本かを知らなきゃだろ？この本の正式なタイトルは「時系列分析と状態空間モデルの基礎」だけど、とりあえず「状態空間モデル」っていうやつの話は本の後半から始まるみたいだから置いといて、「時系列分析」って何かを読み取ろうと思ったんだ。そう思ってまえがきと目次をみると時系列分析とは何か1章に書いてあるらしいから1章を読んだんだけど、「時系列データ」っていうのは具体例から何となくイメージできたんだ。年ごとでも日ごとでも秒ごとでもいいから、一定の時間間隔ごとに1つずつ数字があればいいんだよな？毎月の売上金額とかって…。

ええ、一定の時間間隔ごとに1つずつ数字があればいいですね。もっとも、2つずつでも3つずつでもよいし、数字である必要もないですよ。1秒30フレームの動画があったとして、1フレームごとの静止画像を並べればこれも時系列データといえるでしょう。例えばある曲の1拍ごとのコード進行だって時系列データですよ。静止画像もコードも数字ではありませんけどね。まあ、静止画像もコードもどうにかすれば数値（のベクトル）に対応させることはできますから、時系列データという文脈では既に数字に変換されていることが想定されているとは思いますが。

ジュンの解説逆にややこしいんだけど…。でも、音楽も時系列なのか。なんか親近感湧いてきたな。けど、ここからが問題なんだよ…。時系列データを分析するっていう話になっていくと思うんだけど、最初に出てくる母集団やら母平均の推定っていう話がわからないんだ…。「日本人の身長の平均値」とかいう場合の平均ならわかるけどさ、「2000年1月1日の気温」に母集団も母平均も何もないだろ？「『2000年1月1日』が複数あれば異なる気温だったかもしれない」って言われても、別に複数ないし、あったとしても同じ気温にしかならないかもしれないし…。平均を推定するのが難しいって書いてあるけど、それ以前に平均なんてあるの？

…平均なんてあるのか、ですか。あるかどうかというよりは、あることにしたんです。その前後に「データ生成過程（DGP）」の話は出てきましたか？

あることにした？あ、DGPっていうのは次のページにある。

時系列分析をするときは「時系列データは確率変数列（データ生成過程＝DGP）から生成される」と仮定するんです。つまり、時系列データの各時点の値にはそれを生成する時点ごとの確率分布たちがある、ということにするんです。時系列分析の目的はこの確率分布たち＝DGPの生成過程＝時系列モデルを特定することです。いえ、もう少し正確にいうと、特定した時系列モデルによって「時系列データの将来の値を予測する」とか「複数の時系列を回帰分析する」とか、時系列データからもっと具体的な情報を引き出すことが最終的な目的ですね。逆に、目当ての具体的な情報を得る上では不要な部分を、確率的な揺らぎ＝ノイズとしたのです。本の25ページの言葉を借りれば「未来を予測する情報が含まれていない、純粋な雑音」ですね。ノイズを取り除いて、目的に沿う情報だけを取り出したいんです。確率的なゆらぎがあると決めているから、「『2000年1月1日』が複数あれば異なる気温だったかもしれない」ということになるんです。

ノイズがあると考えて目的の情報を取り出す…？うーん、でも、不要な部分を勝手に考えていいの？ある日の気温が10℃だったというのは、そういう情報じゃん。音楽も時系列って言ってたけど、俺たちの曲に不要な音なんてないだろ？

あくまである目的の上で不要、ですよ。そうですね…あるお店、ドーナツ屋でもいいです、ドーナツ屋が毎週の売上高の成長率をみて、来週の原材料の仕入れを判断しているとしましょう。このお店は売り上げが徐々に伸びていて、成長率の分、原材料を多く仕入れたいんですね。でも、今週はたまたま修学旅行生の集団がこの街にやってきて、ドーナツがとてもよく売れたんです。この週の成長率を参考に来週のために原材料を多く仕入れたら、来週はきっと原材料を余らせてしまうでしょう。あくまで来週の原材料を仕入れる上では、修学旅行生が来たことによる売上増加分はノイズだったんです。音楽でも、例えば4小節などでモチーフを反復する曲は多いでしょう？楽曲のこのような反復の周期を調べるには、4小節ずらして同じ曲を重ねたときに音の波形が重なるかを調べるかもしれません。このとき、モチーフに局所的に入れたアレンジは不要な情報になります。それは周期を調べる上でノイズというだけで、楽曲として不要な音というのではありません。

なるほど…わかってきた気がする。何か調べたいことがあって、それを調べる上で有用な部分と不要な部分に切り分けたいんだな。

もちろん勝手な切り分け方をしてよいわけではありません。有用な部分と不要な部分を切り分ける方法や、切り分け方が適切かどうかを調べる方法を、その本の続きが教えてくれるんでしょう。

でも、何を調べたいのかがまだよくわからないかも。だって、曲に何小節の周期があるかとか聞いて調べればよくない？

そうでしょうか。もしどんな周期性をもつ楽曲が流行するのかをリサーチしたいとして、調査対象が何百曲もあったら1つ1つ聞くのは大変です。楽曲データを入力するたけで判定するモデルがあれば便利だと思います。

あー確かに。

ただそれはあくまで例だったので…僕も詳しくありませんが、音楽データだったら、楽曲へのタグ付け、レコメンデーション、自動作曲などの研究が進んでいるのではないかと思います。「楽曲がどんなジャンルか判定したい」とか、「ある楽曲がある人にオススメできるか判定したい」とか、「勝手にコード進行と主旋律を決定してくれる時系列モデルがほしい」とかでしょうね（これらの目的で必ずしも音楽が時系列として扱われるかはわかりませんが…）。

えー自動作曲！？俺の仕事AIに奪われちゃう！

どうでしょう。確かにコード進行はカデンツのような雛形がありますし、自動作成しやすいと思います。でも、主旋律はもっと自由ですからね。もし本当にハヤトと同じように作曲できるAIをつくりたかったら、そのモデルにはハヤトの作曲の癖が全てプログラミングされていなければならないでしょうけど、そんなことが可能でしょうか。それに、いくらコードとメロディだけ魅力的に生成できても、歌詞やコンセプトがあっての楽曲です。細かなニュアンスを曲に反映するモデルはきっと容易ではないでしょう。作曲を完全に機械の仕事にするのはまだ難しいのではないかと思いますが…。

そっか、よくわかんないけどよかった！

いや僕も知りませんけど…本に戻りましょうか。

うん、時系列分析とは目的にあった時系列モデルを推定することってわかった。2章は時系列データの構造だって。自己相関とコレログラム？

4小節ごとにずっとモチーフを繰り返す楽曲なら、4小節ずらして重ねたときによく重なるということです。そのような曲について、1, 2, 3, 4, 5小節とずらしたときの重なりの度合いを棒グラフにして横に並べていったら、4のところの棒グラフは他より大きくなるでしょう。それがコレログラムです。

なるほど。次のページの外因性っていうのは、まさにさっきのドーナツ屋の修学旅行生のことだな。時系列データの構造って自己相関とかトレンドとか周期性なんだ。28ページの偏自己相関っていうのは自己相関とどう違うんだ…さっきの4小節ごとに繰り返す楽曲なら、8小節ずらしてもよく重なるから、本当に8小節で繰り返す周期性がどれだけあるかを知りたいときは、こっちの偏自己相関をみるってことか…。30ページはホワイトノイズと iid 系列？ホワイトノイズは iid 系列なの？そうじゃないの？

書いてあるように、もし毎時刻の確率分布が正規分布であるようなホワイトノイズであれば、iid 系列です。しかし、一般にはホワイトノイズが独立とは限りません。ホワイトノイズの定義は「毎時刻の期待値はゼロ、毎時刻の分散は一定、異なる時刻の自己相関はゼロ」というだけで独立は課していないんです。毎時刻が正規分布のホワイトノイズの場合は、自己相関がゼロというのが独立と同値になりますが。

じゃあここで iid 系列が出てきたのは何？単なる紹介？「独立性は厳しい仮定です」って言われても、「そっか厳しいんだ」って感じなんだけど…。

そうですね…iid 系列は時系列モデルの最も基本的なパーツとしてつかわれますし、（強）定常という重要な性質をもつので、ここで言及しているのでしょう。31ページで出てくるランダムウォークというモデルが、毎時刻 iid 系列である正規ホワイトノイズを足し上げていくようなモデルになっているでしょう？

iid 系列が重要なパーツになるのか…。でも iid 系列って「毎日サイコロを振った結果」みたいなもんだよな。なんか時系列データっぽくなくない？気温とか売上高とか音楽は、そういうんじゃないじゃん。前の時刻の値がどうだったかといまの値が関係あるよな、たぶん…。

分析対象の時系列データそのものの生成過程が iid 系列ということは少ないかもしれませんが、差分を取ったり、季節変動を除去したりしてデータを定常にしたらその生成過程は iid 系列に近くなるかもしれません。もし対象の売上高データが「『特定の季節変動＋特定のiidノイズ』を足し上げて生成されるもの」だとわかったらうれしいはずです。もう一期先の値の確率分布だってわかります。

差分を取ったりして、定常にする？

第2部 Box-Jenkins 法の内容になりますね。伝統的な時系列分析の手順である Box-Jenkins 法では、まずデータを定常過程に変換することを目指します。定常過程とは何かは35ページから説明がありますが、平均と分散が時刻によらず一定で、自己相関も時点によらず時間差にのみ依存するようなモデルを弱定常過程といいます。つまり、ホワイトノイズは弱定常過程です。自己相関にとどまらず同時分布の形状自体も時間差にのみ依存するようなモデルを（強）定常過程といいます。iid であれば（強）定常です。もっとも、この本では弱定常過程のことを定常過程、そうでない場合を非定常過程と呼んでいるので以降はこれにしたがいましょう。定常過程から生成したデータをプロットすると、ずっと同じ平均値の周りに、ずっと同じばらつきで推移します。他方、非定常過程であったら、平均や分散が拡大していくか、縮小していくか、波打ってしまいます。基本的に時系列モデルは、データを定常過程に変換した上でノイズを説明するように構築します。そのモデルからデータサンプルを生成するときは、ノイズを生成してから逆向きの変換をするんです。

（ノート2があれば）つづく