2020-08-23

Pythonでmode(最頻値)を算出する最良の方法

Python データ分析統計 pandas

Pythonでたまにmode（最頻値）を算出したくなるのですが、いつもどう算出するか迷います。算出方法を整理して、処理時間も計測してみました。

5つの方法
最良の方法（結論）
求め方
処理速度の比較

5つの方法

代表的な方法は以下でしょうか。思いついたものを試してみました。

最良の方法（結論）

どれでも良い。が私が至った結論です。強いて言うならpd.DataFrameを使っていたらそのままpandasを使って、そうでないなら scipyかstatisticsを使う。になると思います。状況に応じて好きなのを使えばいいと思います。

以下の視点で違いを考察しました。

#	視点	差異
1	最頻値が複数ある場合の挙動	あり
2	関数一発か	あり
3	処理速度（最頻値のうち、最小の値を求める場合）	ほぼなし

1. 最頻値が複数ある場合の挙動の差

特に注意しなければ行けないのが、mode（最頻値）が複数ある場合の挙動です。例えば、[1, 2, 2, 3, 3, 4]の最頻値は[2, 3]です。ライブラリや関数によって、以下のように挙動に違いがあるみたいです。statisticsはPythonのバージョンによって挙動が異なるので注意が必要です。

挙動	方法
最頻値を複数返す	pandas, statistics(Python3 $\ge$ 3.8)
最初に出現した値を返す	statistics(Python3 $\ge$ 3.8), collections
最小の値を返す	scipy
例外を返す	statistics(Python3 < 3.8のみ)
自分で処理	numpy

2. 関数一発かの差

numpyとcollectionsは自分で処理を書く必要があります。自由度は高いですが、特別理由がない限りは自分で書く必要もないでしょう。あと最頻値が複数あるケースで、一つに絞り場合は結局maxだかminだか自分で処理する必要はあります。

	方法
関数一発	scipy, pandas, statistics(Python3 $\ge$ 3.8)
自分で処理	numpy, collections

3. 処理速度の差

結果は「scipyと numpyが早い」というのと、驚いたのが「pandasもそれなりに早い」という点です。個人的にはmodeに処理速度を求めることがそんなにないと思いますので、気ににしなくて良いかなと。pd.DataFrameやpd.Sereisを使っている場合はmodeメソッドをそのまま使えばいいと思いますが、modeを算出するためだけにわざわざpd.Seriesとかにキャストするのは個人的に違和感があります。

求め方

メモとして、それぞれの求め方を羅列していきます。

設定

入力する値は以下です。ヒストグラムを見ると、最頻値が[2, 3, 4]であることが分かります。

import matplotlib.pyplot as plt

x = [1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 6]  # has two modes: [2, 3, 4]

# Plot histogram
fig, ax = plt.subplots()
ax.hist(x, bins=np.arange(8) - 0.5, ec='white')
ax.set(xlabel='values', ylabel='Frequency', title='Distribution of x')

f:id:messefor:20200823082040p:plain

1.statistics.mode

標準ライブラリのみで算出することができます。Python3 < 3.8では複数の最頻値があった場合に例外を吐きますので、version3.8+を使うのがおすすめです。

# 最頻値のうち、最初に出現した値を返す@version3.8+
import statistics

statistics.mode(x) # 2

3.8以上では、mode関数は最頻値があった場合、出現した最初のものを返す仕様です。複数最頻値を求めたい場合はmodemulti関数を使いましょう。

If there are multiple modes with the same frequency, returns the first one encountered in the data.

# 複数の最頻値を返したい場合はmodemulti()を使う
statistics.multimode(x) # [2, 3, 4]

例えば、複数の最頻値の中から最小値のみを求めたい場合

# statisticsで複数の最頻値のうち最小の値のみを返す
min(statistics.multimode(x))

2. scipy.stats.mode

scipy.modeは、値と頻度に分けて出力を抽出できます。scipyのmode関数は最頻値があった場合、最小のものを返す仕様です。

If there is more than one such value, only the smallest is returned.

# scipyで最頻値を算出
import scipy
mode, count = scipy.stats.mode(x) # ModeResult(mode=array([2]), count=array([3]))
mode[0] # 2

1つしか最頻値返さないくせにちょっと仰々しい印象を持ちました。

3.numpy.unique

numpyには最頻値を求める関数はなさそうですが、あえて求めるなら以下でしょうか。

# numpyで最頻値を求める
import numpy as np

uniqs, counts = np.unique(x, return_counts=True)
uniqs[counts == np.amax(counts)] # [2, 3, 4]

他のライブラリをインポートする気分じゃないときは良いかもしれません。

複数の最頻値の中から最小値のみを求めたい場合はもちろん以下のようになります

# numpyで複数の最頻値のうち最小の値のみを返す
uniqs[counts == np.amax(counts)].min() # 2

4. pandas.Series.mode

pandasで最頻値を求める場合は以下です。複数の最頻値が返ります。

# pandasで最頻値を求める
import pandas as pd
pd.Series(x, name='values').mode()

# 0    2
# 1    3
# 2    4
# dtype: int64

処理速度を計測すると意外なことに、キャストを含めてもpandas.Series.modeが処理が早かったです。

複数の最頻値の中から最小値のみを求めたい場合は以下のようになります

# pandasで複数の最頻値のうち最小の値のみを返す
pd.Series(x, name='values').mode().min() # 2

5. collections.Counter

最後にcollections.Counterです。複数の最頻値のうち最初の値を返す場合は以下のようになります

from collections import Counter
# collections.Counter複数の最頻値のうち最小の値のみを返す
Counter(x).most_common()[0][0] # 2

ただ複数の最頻値をすべて返したい場合は、あまりCounterを使う利点がないように思いましたので割愛します。

処理速度の比較

処理速度を比較してみました。最頻値を求める際、あまり処理速度をきにすることはないので、個人的には参考程度に思っています。

フェアな比較かどうか微妙なところですが、仕様は以下です。

最頻値の最初のものを一つ求める
入力の1次元配列（リスト）の長さを変えて計測
1条件ごと50回ずつ計測

コードはarticle_mode.py · GitHubこちらにあります。

結果は以下のようになりました。これだけ見るとstatistics, collectionsは入力の長さが大きくなると比較的遅くなるように見えますが、何秒もかかる訳ではないので大差無いと思います。

f:id:messefor:20200823081320p:plain

標準誤差をいれると以下です。

f:id:messefor:20200823081338p:plain

どれでも良い。という結論を出すために長々検証しました。分析としてはコスパが悪くて悪い見本ですが、お役に立てればこれ幸い。日々精進

2020-08-22

Python:venv + pyenvでの環境構築

Python 環境構築仮想環境

一つのOS上で複数バージョンのPythonを切り替えるためのツール、venv +pyenvのメモです。2020年の今、Pythonバージョン3.3以上しか使うことないだろうということで、個人的に venv +pyenvの組み合わせがシンプルで良いかと思います。本投稿では venv +pyenvの役割、導入手順および使い方を整理します。余裕があればpoetryについても触れたいと思います。

前置きが長くなっていますので、venv + pyenvでの仮想環境構築だけ見たい方は導入と使い方、もしくは最後だけ見ればいいかもしれません

実行環境： macOS Mojave

モチベーションと用途の整理
導入と使い方
- venvの導入
  - venvの使い方
- pyenvの導入
  - pyenvで複数のPythonのバージョンをインストール
  - pyenvで複数のPythonのバージョンを切り替え
venv + pyenvで任意のバージョンの仮想環境を構築する
参考文献

モチベーションと用途の整理

まずこれらのツールを導入する際の2つのモチベーションを整理します。たぶん仮想環境という言葉がこの両者を含んでしまっているのが、混乱しやすい原因の一つと思います。

1.パッケージ環境の切り替え

一般にプロジェクトによって必要な外部パッケージ（ライブラリ）は異なります。例えば、画像系ではOpenCVを使うかもしれませんが、画像以外のプロジェクトではなかなか使わないでしょう。不必要なパッケージが多くインストールされている環境は好ましくないため、プロジェクトごとに必要最小限のパッケージがインストールされた環境を用意することになります。このプロジェクトごとのパッケージ環境を切り替えるためのツールの代表がvenvです。

2.Pythonバージョンの切り替え

もうひとつのモチベーションとして、一つのOS上でPythonのバージョンを切り替えて使いたいことがあります。特に以前はPython2系を使っているプロジェクトも多かったため、Python3系とPython2系を切り替えて使う必要がしばしばありました（Googleも2系だったし）。それ以外でも一般に開発を行っていると細かいバージョンの切り替えが必要になります。このPython インタプリタの切り替えを行うツールの代表がpyenvです。

venvとpyenvの用途

使い手から見たvenvとpyenvの主な役割を整理すると、以下のようになります¹。表にもあるように、これらの機能を提供するツールは他にも存在します。

ツール	用途	その他のツール
venv	プロジェクトごとに異なるパッケージをインストールしたいとき	`virtualenv`,`pyenv-virtualenv`, `conda`
pyenv	プロエジェクトごと異なるPythonのバージョンを使いたいとき	`virtualenv`, `conda`

なぜvenv + pyenvか

どのツールを選ぶかは、主に次の要素に依存する気がします。

使いたいPythonのバージョン
使っているOS
好み

Python3.3以上では、パッケージ環境（仮想環境）の切り替えとしてvenvが標準で使えます。Pythonのバージョンを一つしか使ってない方はこのvenvを使えば十分です。たとえばPythonはOSにインストールされた3.6を使って、プロジェクトによってパッケージ群だけ切り替えるケースがこれに当たります。開発しているわけでは無ければ、ほぼこれで事足りると思います。同様に複数のPythonバージョンを切り替えたくなった場合でも、Python3.3以上のバージョンの切り替えであれば、venv に加えてPython インタプリタ切り替えツールとして pyenvを使えば良いと思います。私はこのケースに該当します。

一方でバージョン3.3未満のPythonを使う場合はvenvが使えないので、pyenv+ pyenv-virtualenvを使うのが良いでしょう。pyenv-virtualenvはpyenvのプラグインで、仮想環境の切り替え機能を提供します。

結局、ニーズに合わせて下表のようにまとまる気がします。

	Python3.3以上のみ使う	Python3.3未満も使う
パッケージ環境のみ切り替えたい	`venv`	`virtualenv`（？しらん）
パッケージ環境とPythonバージョンを切替えたい	`pyenv`+`venv`	`pyenv`+`pyenv-virtualenv`

一昔前は、pyenvがWindowsに対応していなかったことと、データ分析界隈でAnacondaがよく使われていたこともあり、Windowsではcondaによるパッケージやインタプリタ切り替えを行っていました。現在はpyenv-winというのがあるらしく状況が変わっているようです。私も以前WindowsではAnacondaを使っていましたが、環境が管理しにくい印象があって現在は使わなくなりました（というかWindowを使っていない）

導入と使い方

venv + pyenvで十分である。という勝手な結論を前提にして、それらの導入方法と使い方を簡単に説明します。

venvの導入

Python3.3以上ではvenvは標準ライブラリに含まれているので、導入の必要はありません。

venvの使い方

仮想環境の一般的な利用手順は次のようになります

プロジェクトごとに仮想環境を作成
仮想環境に入る
環境構築や開発を行う
1. ライブラリのインストール、環境変数設定など
2. コーディング、実行
仮想環境を抜ける

1.プロジェクトごとに仮想環境を作成

作成する前に、念のためシステムにインストールされているPythonのバージョンを確認しましょう。

# pythonのバージョン確認
$ python -V

>>> Python 3.8.5

ここでバージョンが3.3以上ならば問題ないのですが、それ未満だとvenvは使えないので3.3以上をインストールする必要があります。インストール方法は"python3 インストール"とかでググってみてください。

バージョン3.3以上でしたら、pythonからvenvモジュールを実行して仮想環境を作成できます。このとき同時に仮想環境の情報を保存するためのディレクトリパスを指定します。

# .venv/というディレクトリに新しい仮想環境を構築
$ python -m venv .venv

ここでは.venvというディレクトリ名を指定しているので、上記コマンドを実行すると.venv/がカレントディレクトリに作成されます。

2.仮想環境に入る

上で指定した.venv/内にあるactivateコマンドを使って仮想環境に入ります。仮想環境に入るとプロンプトの前に(.venv)という表示が現れるので、環境に入ったことが分かります。

# 仮想環境に入る
$ source .venv/bin/activate

# (.env)という表示がプロンプトの前に現れる
(.venv)$

現在の環境は作成したばかりのまっさらな仮想環境なので、外部パッケージは何もインストールされていません。ためにしnumpyを呼び出してみるとエラーが出ることが分かります。

# pythonのインタラクティブモードに入る
(.venv)$ python

Python 3.8.5 (default, Aug 22 2020, 07:29:50)
[Clang 10.0.1 (clang-1001.0.46.4)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>

# numpy をインポートしてみるとエラーが出る
>>> import numpy as np  
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'numpy'

3.環境構築や開発を行う

このまっさらな環境に必要なパッケージをインストールし、コードを書いて、実行していきます。

パッケージのインストール

パッケージのインストールはpipでできます。仮想環境に入った状態でインストールしたパッケージは仮想環境内のみで利用可能です。外の環境に影響を及ぼしません。

# numpyとjupyterをインストール
(.venv)$pip install numpy jupyter

Pythonの実行

仮想環境に入った状態で、pythonを実行すれば良いです。ファイルを指定したりして、普通に実行してください。

# いつもどおりpythonを実行すれば良い
(.venv)$ python example.py

Jupyter notebookを使う

仮想環境に入った状態でjupyter notebookを立ち上げれば、仮想環境下でnotebookを走らせることができます。

(.venv)$ jupyter notebook

Hydrogenで使う（kernelspecをインストール）

仮想環境に入った状態でkernelspecをインストールすれば、AtomのHydrogenや他のJupyter notebookからこの仮想環境を選択できるようになります。

(.venv)$python -m ipykernel install --user --name kernel_name

4.仮想環境を抜ける

deactivateコマンドで仮想環境から抜けられます。deactivateとだけ打てばOKです。プロンプトの前の表示こが普段どおりに戻れば、環境から抜けられています。

# 仮想環境に入る
(.venv)$ deactivate

# (.env)という表示がプロンプトが消える
$

pyenvの導入

つぎはPythonバージョンの切り替えツールpyenvです。pyenvは環境変数PATHを必要に応じて変更することでPython インタプリタを意図通りに切り替えてくれます。pyenvのインストールはpyenvのgithubページの手順に沿って行えば簡単にできます。

ここではGithub経由、かつbashを使っている環境でのでインストール例を本家ページから引用します。環境が異なる方は本家ページを見てください。

# pyenvのインストール。clone先は$HOME/.pyenvがおすすめ
git clone https://github.com/pyenv/pyenv.git ~/.pyenv

# 環境変数の設定
echo 'export PYENV_ROOT="$HOME/.pyenv"' >> ~/.bash_profile
echo 'export PATH="$PYENV_ROOT/bin:$PATH"' >> ~/.bash_profile

# PATH挿入と自動補完の設定。 eval "$(pyenv init -)"がきちんと.bash_profileに記載されたか確認しよう
echo -e 'if command -v pyenv 1>/dev/null 2>&1; then\n  eval "$(pyenv init -)"\nfi' >> ~/.bash_profile

# シェルの再起動
exec "$SHELL"

pyenvで複数のPythonのバージョンをインストール

pyenvを使って任意のPythonバージョンをインストールしてみましょう。pyenv install --listで利用可能なPythonのバージョン一覧を確認できます。anacondaやminicondaもインストールできるのが面白いですね。

# pyenvで利用可能なPythonバージョン一覧を確認
$ pyenv install --list

Available versions:
  2.1.3
  2.2.3
  2.3.7
  2.4.0
...
  3.8.3
  3.8.4
  3.8.5
  3.9.0b5
  3.9-dev
  3.10-dev
  activepython-2.7.14
  activepython-3.5.4
  activepython-3.6.0
  anaconda-1.4.0
  anaconda-1.5.0
  anaconda-1.5.1
...

ここでは3.8.3をインストールしてみたいと思います。インストールにはちょっと時間がかかります。

# Python3.8.3をインストール
$ pyenv install 3.8.3

インストール完了後、pyenv versionsとすれば、現在の環境にインストールされているPython インタプリタ一覧が確認できます。

# インストールされているPython一覧を確認する
$ pyenv versions
* system (set by /Users/daisuke/.pyenv/version)
  3.8.3

上の例では、system と記載されているものと3.8.3がインストールされているのが確認できます。 * マークがついているのが現在設定されているPythonのバージョンです。systemというのは元からシステムにインストールされているpythonコマンドを表します。

pyenvで複数のPythonのバージョンを切り替え

pyenvでのバージョンの切り替え方法は、いくつかあるのですが、ここでは現在のシェルで一時的にバージョンを切り替える方法pyenv shellを説明します。

特定のバージョン（例えば、3.8.3）に切り替えたい場合は、pyenv shell 3.8.3で可能です。元のバージョンに戻りたいときは、shellを終了するか、同じコマンドでsystemを指定すれば良いかと思います。

# バージョン3.8.3に切り替える
$ pyenv shell 3.8.3

＃ 現在のバージョンを確認する
$ pyenv version
3.8.3 (set by PYENV_VERSION environment variable)

# バージョン3.8.3から元にもどす
$ pyenv shell system

venv + pyenvで任意のバージョンの仮想環境を構築する

今回これを書きたかっただけなのですが、前置きが長かった。。。。

venv+pyenvを使って、任意のバージョン仮想環境（パッケージ環境）を構築したいときは以下のようにします。単純に構築したいPythonのバージョンに切り替えて、仮想環境を作成するだけです。

# pyenvでPythonバージョンを切り替える
$ pyenv shell 3.8.3

# venvで新しい仮想環境を構築
$ python -m venv .venv

これにより、指定のPythonバージョンの仮想環境が.venv/に構築されます。仮想環境に入って実際に確認すると意図通りのPythonバージョンに設定されていることが分かります。

# 仮想環境に入る
$ source .venv/bin/activate

# 仮想環境内のPythonのバージョン確認
(.venv)$ python -V
Python 3.8.3

ちなみに仮想環境内のPythonとpyenv shellした直後のPythonは同じバージョン（3.8.3）ですが、2つのコマンドの実態は別物になります。

下記のように仮想環境内では.venv/の中にインストールされたPython インタプリタを参照しています。

# 仮想環境のPythonは.venv/内のpythonを見にいっている
(.evn)$ which python
/Users/messefor/prj/blog200822/.venv/bin/python

一方でpyenv shellによる切り替え直後はpyenvによってPATHに挿入されたPython インタプリタを参照しています。shimsというのがpyenvによって挿入された参照先になります

$ pyenv shell 3.8.3

$ which python
/Users/messefor/.pyenv/shims/python

久々にPythonの環境構築の話題を書いたら長くなりました。。。。次はpoetryについてまとめたいと思います。日々精進

参考文献

実際にはvenvやvirtualenvが提供するのは仮想環境切り替えで、Python インタプリタなども含めた環境を指しているかと思いますが、ここではパッケージの切り替え、Pythonバージョンの切り替えの二つで分類しています。↩

2020-08-14

時系列分析：ARモデルの性質と定常性（1/2）

データ分析時系列分析確率過程 Python

ARモデルはシンプルな時系列モデルですが、自己共分散や自己相関はMAモデルほど簡素な表現はなく、Yule-Walker方程式のような連立方程式を解いていくのが一般的のようです。本投稿では2回に分けてARモデルの期待値、自己相関および定常性などの性質を整理したいと思います。

ARモデル
- AR(1)の例
期待値
- AR(1)の場合：
- AR(p)の場合：
自己分散・共分散
- AR(1)の場合：
- AR(2)の場合：
  - 分散：
  - 自己共分散：
まとめ
参考文献

ARモデル

ARモデルは現在の値が、過去の自分自身の状態を引きずっている（記憶している）という構造を持つモデルです。世の中の現象ってだいたい過去の状態を引きずっていそうです。 ARモデルではこの過去の記憶を、過去の値の線形和で表現します。p次のARモデルは以下のように定式化されます。

$\displaystyle{ \begin{align} y_t &= c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + \dots + \phi_{p}y_{t-p} + \varepsilon_t\\ &=c + \sum_{k=1}^{p}\phi_ky_{t-k} + \varepsilon_t\\ \\ \\ \varepsilon_t &\sim \operatorname{W.N.}(\sigma_{\varepsilon}^2) \end{align} }$

ここで $\displaystyle{c}$ は定数項でベースのレベルを決めます（直流成分）。中心化をして $\displaystyle{c=0}$ を前提とすることも多いです。 $\displaystyle{\varepsilon_t}$ はホワイトノイズで、 $\displaystyle{\varepsilon_t}$ は、 $\displaystyle{t-1}$ 以前には含まれていない $\displaystyle{t}$ 時点で新しく加わった情報でイノベーションとも呼ばれています。

AR(1)の例

いくつかAR過程の具体例を見てみます。ここでは $\displaystyle{c=0}$ とし、 $\displaystyle{\phi_1}$ の値をいくつか変えて生成し、その原系列とACFをプロットします。生成に使ったPython スクリプトはw3_1.py · GitHubに置いています。

$\displaystyle{\phi_1=1}$ の場合(ランダムウォーク ):

1次のAR過程のうち $\displaystyle{\phi_1=1}$ かつホワイトノイズに同分布を仮定するとランダムウォークと呼ばれる有名な過程になります。

$\displaystyle{ y_{t} = y_{t-1} + \varepsilon_t \\ \\ \varepsilon_t \sim \operatorname{iid}(0, \sigma_{\varepsilon}^2) }$

生成したランダムウォークの推移とACFは次のようになりました。

f:id:messefor:20200814182540p:plain — acf_phi_0.3

$\displaystyle{\phi_1=0.3}$ の場合:

f:id:messefor:20200814182500p:plain

$\displaystyle{\phi_1=1.1}$ の場合:

f:id:messefor:20200814182618p:plain

同じAR(1)過程でも $\displaystyle{\phi}$ のとる値によって印象は随分違います。 $\displaystyle{\phi_1=1.1}$ の例が発散していることからも明らかなようにAR過程は $\displaystyle{\phi}$ の値によって定常性も異なります。

期待値

AR過程の期待値を計算してみます。

AR(1)の場合：

1次のARモデル $\displaystyle{y_t = c + \phi_{1}y_{t-1} + \varepsilon_t}$ の期待値をとると

$\displaystyle{ \begin{align} E[y_{t}] &= E[c + \phi_{1}y_{t-1} + \varepsilon_t]\\ &=c + \phi_{1}E[y_{t-1}] \end{align} }$

ここで定常を仮定とすると期待値は時刻 $\displaystyle{t}$ に依存しないため、 $\displaystyle{E[y_t]=E[y_{t-1}]}$ となります。 $\displaystyle{E[y_t]=E[y_{t-1}]=\mu}$ と置くと次のようになります。

$\displaystyle{ \begin{align} \mu &=c + \phi_{1}\mu \\ \\ \mu &= \frac{c}{1-\phi_1} \end{align} }$

AR(p)の場合：

p次の場合も同様に期待値をとって

$\displaystyle{ \begin{align} E[y_t] &= E[c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + \dots + \phi_{p}y_{t-p} + \varepsilon_t]\\ \\ &=c + \phi_{1}E[y_{t-1}] + \phi_{2}E[y_{t-2}] + \dots + \phi_{p}E[y_{t-p}] \\ \end{align} }$

定常を仮定して、 $\displaystyle{E[y_t]=E[y_{t-1}]=\dots=\mu}$ と置くと次のようになります。

$\displaystyle{ \begin{align} \mu &=c + \phi_{1}\mu + \phi_{2}\mu +\dots + \phi_{p}\mu \\ \\ \mu &= \frac{c}{1-\phi_{1}-\phi_{2} - \dots - \phi_{p}} \end{align} }$

また上式より直流成分 $\displaystyle{c=0}$ のとき、定常なARモデルの期待値は常に0になることが分かります。

自己分散・共分散

自己分散・共分散は少し工夫が必要になります。まずは愚直に計算してみましょう。

AR(1)の場合：

分散：

AR(1)モデル $\displaystyle{y_t = c + \phi_{1}y_{t-1} + \varepsilon_t}$ の分散を計算します。

$\displaystyle{ \begin{align} \operatorname{Var}[y_{t}] &= \operatorname{Var}[c + \phi_{1}y_{t-1} + \varepsilon_t]\\ &= \operatorname{Var}[c] + \operatorname{Var}[\phi_1y_{t-1}] + \operatorname{Var}[\varepsilon_t] \\&+ 2\operatorname{Cov}(c, \phi_1y_{t-1})+ 2\operatorname{Cov}(\phi_1y_{t-1}, \varepsilon_t) + 2\operatorname{Cov}(c, \varepsilon_t) \\ &=\phi_{1}^2\operatorname{Var}[y_{t-1}] + \sigma_{\varepsilon}^2 \end{align} }$

定常を仮定して、 $\displaystyle{\operatorname{Var}[y_t]=\operatorname{Var}[y_{t-1}]=\gamma_0}$ と置くと以下のようになります。

$\displaystyle{ \begin{align} \gamma_0 &=\phi_{1}^2\gamma_0 + \sigma_{\varepsilon}^2 \\ \\ \gamma_0 &= \frac{\sigma_{\varepsilon}^2}{1-\phi_1^2} \tag{1} \end{align}\\ }$

自己共分散：

$\displaystyle{y_t}$ と $\displaystyle{y_{t-k}}$ の共分散をとります。

$\displaystyle{ \begin{align} \operatorname{Cov}(y_{t}, y_{t-k}) &= \operatorname{Cov}(c + \phi_{1}y_{t-1} + \varepsilon_t, y_{t-k})\\ &= \operatorname{Cov}(c, y_{t-k}) + \phi_{1}\operatorname{Cov}(y_{t-1}, y_{t-k}) + \operatorname{Cov}(\varepsilon_t, y_{t-k})\\ &= \phi_{1}\operatorname{Cov}(y_{t-1}, y_{t-k}) \end{align} }$

2行目への展開は、共分散を分配しています。3行目への展開は、定数や独立な変数との共分散は0であることを利用しています。分散・共分散の計算は以下を参考にしてください。

messefor.hatenablog.com

ここで $\displaystyle{\operatorname{Cov}(y_{t}, y_{t-k})}$ はラグ $\displaystyle{k}$ の自己共分散なので $\displaystyle{\gamma_k}$ 、また $\displaystyle{y_{t-1}}$ と $\displaystyle{ y_{t-k}}$ のラグは $\displaystyle{k-1}$ なので同様に $\displaystyle{\operatorname{Cov}(y_{t-1}, y_{t-k})}$ を $\displaystyle{\gamma_{k-1}}$ と置くと

$\displaystyle{ \gamma_k = \phi_1\gamma_{k-1} \tag{2} }$

という自己共分散の関係式が現れます。 $\displaystyle{k=1}$ の自己共分散は $\displaystyle{(1)(2)}$ 式を使って次のように表せます。

$\displaystyle{ \begin{align} \gamma_1 = \phi_1\gamma_0 &=\phi_1\frac{\sigma_{\varepsilon}^2}{1-\phi_1^2} \end{align} }$

Yule-Walker方程式

さらに $\displaystyle{(2)}$ 式を $\displaystyle{\gamma_0}$ で割ると、自己相関でも同様の関係が成り立つことが分かります。

$\displaystyle{ \begin{eqnarray} \frac{\gamma_k}{\gamma_0} = \phi_1\frac{\gamma_{k-1}}{\gamma_0} \\ \rho_k = \phi_1\rho_{k-1} \tag{3} \end{eqnarray} }$

$\displaystyle{(3)}$ 式はYule-Walker方程式と呼ばれているものです。

AR過程の自己相関が $\displaystyle{y_t}$ が従うAR過程と同一の係数をもつ差分方程式に従うことを示す

引用：経済・ファイナンスデータの計量時系列分析- 沖本竜義

$\displaystyle{(3)}$ 式の関係性を繰り返し適用すると、ラグ $\displaystyle{k}$ の自己相関は $\displaystyle{\rho_0}$ を使って次のように書けます。

$\displaystyle{ \rho_k = \phi_1\rho_{k-1} = \phi_1^2\rho_{k-2} = \dots = \phi_1^k\rho_0 }$

ラグ0の自己相関 $\displaystyle{\rho_0=1}$ なので

$\displaystyle{ \rho_k = \phi_1^k }$

このようにAR(1)の自己相関は非常にシンプルな形で表現できることが分かります。また $\displaystyle{|\phi_1| \lt 1}$ のとき自己相関は $\displaystyle{k}$ が大きくなるにつれて減衰することが分かります。

AR(2)の場合：

AR(1)のケースでは素直に分散・自己共分散を算出できましたが、同じことをAR(2)でも試してみます。

分散：

AR(2)モデル $\displaystyle{y_t = c + \phi_{1}y_{t-1}+ \phi_{2}y_{t-2} + \varepsilon_t}$ の分散を計算します。

$\displaystyle{ \begin{align} \operatorname{Var}[y_{t}] &= \operatorname{Var}[c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2}+ \varepsilon_t]\\ &= \operatorname{Var}[\phi_1y_{t-1}] +\operatorname{Var}[\phi_2y_{t-2}]+\operatorname{Var}[\varepsilon_t] + 2\operatorname{Cov}(\phi_1y_{t-1}, \phi_2y_{t-2}) \\ &=\phi_{1}^2\operatorname{Var}[y_{t-1}] + \phi_{2}^2\operatorname{Var}[y_{t-2}] + \sigma_{\varepsilon}^2 + 2\phi_{1}\phi_{2}\operatorname{Cov}(y_{t-1}, y_{t-2}) \end{align} }$

定常を仮定すると $\displaystyle{\operatorname{Var}[y_{t}] = \operatorname{Var}[y_{t-1}] = \operatorname{Var}[y_{t-1}] = \gamma_0, \operatorname{Cov}(y_{t-1}, y_{t-2})=\gamma_1}$ といえるので

$\displaystyle{ \gamma_0 =\phi_{1}^2\gamma_0 + \phi_{2}^2\gamma_0 + \sigma_{\varepsilon}^2 + 2\phi_{1}\phi_{2}\gamma_1 \tag{4} }$

このように $\displaystyle{\gamma_0}$ と $\displaystyle{\gamma_1}$ の関係式が出てきました。つまり分散を知るには、ラグ1の自己共分散 $\displaystyle{\gamma_1}$ が必要になります。

自己共分散：

ラグ1の自己共分散 $\displaystyle{\gamma_1}$ を計算します。

$\displaystyle{ \begin{align} \gamma_1 &= \operatorname{Cov}(y_{t}, y_{t-1})\\ &= \operatorname{Cov}(c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2}+ \varepsilon_t, y_{t-1})\\ &=\operatorname{Cov}(c, y_{t-1}) + \phi_{1}\operatorname{Cov}(y_{t-1},y_{t-1}) + \phi_{2}\operatorname{Cov}(y_{t-2}, y_{t-1})+ \operatorname{Cov}(\varepsilon_t, y_{t-1})\\ &=\phi_{1}\gamma_0 + \phi_{2}\gamma_1 \\ \\ \gamma_1 &= \frac{\phi_1}{1-\phi_2}\gamma_0 \tag{5} \end{align} }$

$\displaystyle{(4)}$ 式に $\displaystyle{(5)}$ を代入して整理すると、

$\displaystyle{ \begin{eqnarray} \frac{1-\phi_2 - \phi_2^2 + \phi_2^3 - \phi_1^2 - \phi_1^2\phi_2}{1-\phi_2}\gamma_0 = \sigma_{\varepsilon}^2 \\ \gamma_0 = \frac{1-\phi_2}{1+\phi_1}\frac{\sigma_{\varepsilon}^2}{((1-\phi_2)^2 - \phi_1^2)} \tag{6} \end{eqnarray} }$

ということで分散が求まりました。さらに $\displaystyle{\gamma_1}$ は $\displaystyle{(5)}$ 式に $\displaystyle{(6)}$ を代入して

$\displaystyle{ \gamma_1 = \frac{\phi_1}{1+\phi_1}\frac{\sigma_{\varepsilon}^2}{((1-\phi_2)^2 - \phi_1^2)} }$

となります。同様に $\displaystyle{\gamma_2}$ 以降も求めることができます。このようにAR過程ではあるラグでの自己共分散を求める際、他のラグの自己共分散の値を使って連立方程式を解く必要があります。これらの操作は次のYule-Walker方程式を使って統一的かつ簡潔に表現できます。

$\displaystyle{ \rho_k = \phi_1\rho_{k-1} + \phi_2\rho_{k-2} + \dots + \phi_p\rho_{k-p}, \ k >0 }$

次回、Yule-Walker方程式を用いたAR(p)の自己相関関数や定常性について整理します。

まとめ

ARモデルは過去の状態の線形和使って現在の値を表現します
ARモデルの定常性は係数によって決まります
AR(1)の期待値、分散、自己相関係数はそれぞれ次のようになります

$\displaystyle{ \begin{align} \mu &= \frac{c}{1-\phi_1}\\ \\ \gamma_0 &= \frac{\sigma_{\varepsilon}^2}{1-\phi_1^2}\\ \\ \rho_k &= \phi_1^k \end{align}\\ }$

AR過程の自己相関は $\displaystyle{y_t}$ が従うAR過程と同一の係数をもつ差分方程式に従います（Yule-Walker方程式）

$\displaystyle{ \rho_k = \phi_1\rho_{k-1} }$

ARモデルの分散や自己相関は上記の関係性を使うことで求めることができます

AR(2)の共分散の式展開が合わずに時間を無駄に使いました。とほほ。中年の残りすくない時間を。。。日々精進

参考文献

自己回帰モデル - Wikipedia
time series - Variance of a stationary AR(2) model - Cross Validated
Time Series AR(2) Variance Derivation : AskStatistics
経済・ファイナンスデータの計量時系列分析- 沖本竜義

2020-08-12

時系列分析：MAモデルの反転可能性

ARMAモデルデータ分析時系列分析確率過程

自己相関構造は時系列モデル選択の重要な基準になります。しかしながらMAモデルでは同じ自己相関構造を持つモデル（パラメタ）が複数存在するため、その中からMAモデルを一つ選ぶ基準が必要です。この基準として反転可能性の議論がでてきます。本投稿ではMAモデルの反転可能性について確認します。時系列分析において反転可能性は重要な性質の一つで、ARモデルの収束性判定でも用いられます。

MAモデルの選択
反転可能性（invertiblity）
まとめ

MAモデルの選択

自己相関構造は時系列モデル選択の目安となりますが、MAモデルでは同じ期待値や自己相関構造を持つパラメタは複数存在します。以下の２つのMA(1)モデルそれぞれの自己相関構造を確認してみましょう。

$\displaystyle{ \begin{eqnarray} y_t = \varepsilon_t + 2\varepsilon_{t-1} \tag{1}\\ y_t = \varepsilon_t + \frac{1}{2}\varepsilon_{t-1} \tag{2}\\ \\ \varepsilon_t \sim \operatorname{W.N}(\sigma^2) \end{eqnarray} }$

$\displaystyle{(1)}$ の期待値と自己共分散を定義どおりに計算すると

$\displaystyle{ \begin{align} \mu &= E[y_t] = E[\varepsilon_t + 2\varepsilon_{t-1}] = 0 \\ \\ \gamma_k &= \operatorname{Cov}[\varepsilon_t + 2\varepsilon_{t-1}, \varepsilon_{t+k} + 2\varepsilon_{t+k-1}]\\ &= E[(\varepsilon_t + 2\varepsilon_{t-1})(\varepsilon_{t+k} + 2\varepsilon_{t+k-1})] - E[\varepsilon_t + 2\varepsilon_{t-1}]E[\varepsilon_{t+k} + 2\varepsilon_{t+k-1}] \\ &= E[\varepsilon_t\varepsilon_{t+k}] + 2E[\varepsilon_t\varepsilon_{t+k-1}] + 2E[\varepsilon_{t-1}\varepsilon_{t+k}] + 4E[\varepsilon_{t-1}\varepsilon_{t+k-1}] \tag{3} \end{align} }$

ホワイトノイズの性質から時点が一致する自己共分散以外は0になるので、 $\displaystyle{(3)}$ 式はラグ $\displaystyle{k}$ によって場合分けされ

$\displaystyle{ \gamma_k = \left\{ \begin{array}{ll} 5\sigma^2 & k = 0\\ 2\sigma^2 & k = 1\\ 0 & k > 1\\ \gamma(-k) & k < 0 \end{array} \right. }$

　これより自己相関は

$\displaystyle{ \rho_k = \frac{\gamma_k}{\gamma_0} = \left\{ \begin{array}{ll} 1 & k = 0\\ \frac{2}{5} & k = 1\\ 0 & k > 1\\ \rho(-k) & k < 0 \end{array} \right. }$

なお一般にMA(q)モデルの自己相関は次のようになります。これが分かっていれば、上記はすぐに算出できます。

$\displaystyle{ \rho_k = \frac{\sum_{i=0}^{q-k}\beta_i\beta_{i+k}}{\sum_{i=0}^{q}\beta_i^2} }$

MA(q)モデルの自己相関の考え方については以下を参考にしてください。

messefor.hatenablog.com

同様に $\displaystyle{(2)}$ の期待値、自己相関を計算すると

$\displaystyle{ \begin{align} \mu &= E[y_t] = E[\varepsilon_t + \frac{1}{2}\varepsilon_{t-1}] = 0 \\ \\ \gamma_k &= \operatorname{Cov}[\varepsilon_t + \frac{1}{2}\varepsilon_{t-1}, \varepsilon_{t+k} + \frac{1}{2}\varepsilon_{t+k-1}]\\ &= E[(\varepsilon_t + \frac{1}{2}\varepsilon_{t-1})(\varepsilon_{t+k} + \frac{1}{2}\varepsilon_{t+k-1})] - E[\varepsilon_t + \frac{1}{2}\varepsilon_{t-1}]E[\varepsilon_{t+k} + \frac{1}{2}\varepsilon_{t+k-1}] \\ &= E[\varepsilon_t\varepsilon_{t+k}] + \frac{1}{2}E[\varepsilon_t\varepsilon_{t+k-1}] + \frac{1}{2}E[\varepsilon_{t-1}\varepsilon_{t+k}] + \frac{1}{4}E[\varepsilon_{t-1}\varepsilon_{t+k-1}] \end{align} }$

$\displaystyle{ \gamma_k = \left\{ \begin{array}{ll} \frac{5}{4}\sigma^2 & k = 0\\ \frac{1}{2}\sigma^2 & k = 1\\ 0 & k > 1\\ \gamma(-k) & k < 0 \end{array} \right. \\ }$

$\displaystyle{ \\ \rho_k = \frac{\rho_k}{\rho_0}\left\{ \begin{array}{ll} 1 & k = 0\\ \frac{2}{5} & k = 1\\ 0 & k > 1\\ \rho(-k) & k < 0 \end{array} \right. }$

２つのモデルの期待値、自己相関は全く同じになります。沖本本によると

同一の期待値と自己相関構造をもつMA(q)過程は一般的に $\displaystyle{2^ q}$ 個あることが知られている

とのことなので、自己相関構造からMAモデルを一つ選ぶ基準が必要です。

反転可能性（invertiblity）

この基準として登場するのが反転可能性（invertiblity)です。MAモデルは条件によって、以下のようにAR( $\displaystyle{\infty}$ )で表すことができます。これを反転と呼んでいます。

$\displaystyle{ \begin{eqnarray} y_t = \sum_{k=1}^{\infty}\pi_ky_{t-k} + \varepsilon_t \\ \mathrm{or} \\ \varepsilon_t = y_t - \sum_{k=1}^{\infty}\pi_{k}y_{t-k} \tag{4} \end{eqnarray} }$

$\displaystyle{(4)}$ の形を見ると分かりやすいですが、AR( $\displaystyle{\infty}$ )に反転させた形でみると $\displaystyle{\varepsilon_t}$ が $\displaystyle{y_t}$ と、過去の $\displaystyle{y}$ の線形和の差になっています。これは現在の $\displaystyle{y_t}$ のうち、過去の $\displaystyle{y}$ を使っても表現できない部分が $\displaystyle{\varepsilon_t}$ である、と解釈できます。反転可能なMA過程は、解釈が容易になるなど良い性質を持っています。

MAモデルの反転のイメージ

なぜMAモデルがAR( $\displaystyle{\infty}$ )で表せるのでしょうか。MA(1)モデルの例で確認します。ここでは等比数列の部分和を用いて導出してみます。まずMA(1)モデルをラグオペレータ（後退オペレータ） $\displaystyle{B}$ の多項式を用いて表します。

$\displaystyle{ \begin{eqnarray} y_t &= \varepsilon_t + \beta_1\varepsilon_{t-1}\\ &=\varepsilon_t - \beta_1B\varepsilon_t \\ &=(1 - \beta_1B)\varepsilon_t\\ &=\beta(B)\varepsilon_t\\ \mathrm{where}\\ \\ \beta(B) &= 1 - \beta_1B \tag{5} \end{eqnarray} }$

ラグオペレータについては以下の記事を参考にしてください。

messefor.hatenablog.com

ここで $\displaystyle{\beta(B)}$ はラグオペレータの多項式です。ここでもし $\displaystyle{\beta(B)}$ の逆演算を定義できるなら次のように撹乱項を $\displaystyle{y_t}$ とラグオペレータで表現できます。

$\displaystyle{ \begin{eqnarray} \varepsilon_t &= \beta(B)^{-1}y_t \tag{6} \end{eqnarray} }$

実際に $\displaystyle{B}$ は単位円上の複素数と考えて良いようで、 $\displaystyle{\beta(B)^ {-1}}$ は次のようになります。

$\displaystyle{ \begin{eqnarray} \beta(B)^{-1} &= \frac{1}{1 - \beta_1B} \tag{7} \end{eqnarray} }$

ここで高校数学の等比数列の部分和から以下の関係性を使います。

$\displaystyle{ \begin{align} S_n &= a + ar + ar^2 + \dots + ar^{n-1}\\ &= a\frac{1-r^n}{1-r} \tag{8}\\ \end{align} }$

$\displaystyle{n}$ を $\displaystyle{\infty}$ にとばした場合、 $\displaystyle{(8)}$ 式は $\displaystyle{|r| \lt 1}$ のとき収束し次のようになります

$\displaystyle{ \begin{align} \lim_{n \to \infty}S_n &= a + ar + ar^2 + \dots + ar^{n-1} + \dots\\ &= a\frac{1}{1-r} \tag{9}\\ \end{align} }$

この関係性を使って $\displaystyle{(7)}$ 式を展開してみます、 $\displaystyle{(8)}$ 式で $\displaystyle{a=1}$ 、 $\displaystyle{r=\beta_1B}$ と考えて、 $\displaystyle{|r| \lt 1}$ のとき

$\displaystyle{ \begin{align} \beta(B)^{-1} &= \frac{1}{1 - \beta_1B} \\ &=1+\beta_1B + \beta_1^2B^2 + \dots + \beta_1^{n-1}B^{n-1} + \dots\\ \end{align} }$

さらにこれを用いて、 $\displaystyle{(6)}$ 式を展開すると

$\displaystyle{ \begin{align} \varepsilon_t &= \beta(B)^{-1}y_t\\ &=(1+\beta_1B + \beta_1^2B^2 + \dots + \beta_1^{n-1}B^{n-1} + \dots)y_t\\ &=y_t + \beta_1y_{t-1} + \beta_1^2y_{t-2} + \dots + \beta_1^{n-1}y_{t-n+1} + \dots\\ &=y_t + \sum_{k=0}^{\infty}\beta_1^ky_{t-k} \end{align} }$

が成り立ちます。これで $\displaystyle{(4)}$ 式と同じ形になりました。

ここで $\displaystyle{|r| \lt 1}$ の条件は $\displaystyle{|r|=|\beta_1B| \lt 1}$ と、 $\displaystyle{B}$ が単位円上の複素数であることを考えると $\displaystyle{|B| = 1}$ のため、最終的に $\displaystyle{|\beta_1| \lt 1}$ が満たすべき条件となります。

MA過程の反転可能条件

MA(1)の反転を行う過程では、 $\displaystyle{|r| \lt 1}$ という条件が必要ということがわかりました。一般にMA(q)過程の反転可能性の条件はMA特性方程式の根を使って設定されます。MA特性方程式とは、ラグオペレータ多項式のラグオペレータ $\displaystyle{B}$ を複素数 $\displaystyle{z}$ と置き換え、多項式を0とした方程式で次の形をとります。

$\displaystyle{ \begin{align} \theta(z) &= 1 + \theta_1z + \theta_2z^2 + \dots + \theta_qz^q = 0 \end{align} }$

MA過程の判定可能性条件:

MA特性方程式 $\displaystyle{\theta(z)}$ の根がすべて単位円の外に存在する（ $\displaystyle{z}$ の解の絶対値が1より大きい）場合、そのMA過程は判定可能

この反転可能性の制約を入れることで、期待値と自己相関構造からMAモデルは一つに定まります。

MA(1)の反転可能性

先程のMA(1)の例で反転可能性条件を確認してみます。 MA(1)モデル $\displaystyle{y_t = \varepsilon_t + \beta_1\varepsilon_{t-1}}$ のMA特性方程式は、 $\displaystyle{(5)}$ 式のラグオペレータ多項式 $\displaystyle{\beta(B) = 1 - \beta_1B}$ を考えて

$\displaystyle{ \begin{eqnarray} \beta(z) = 1 - \beta_1z = 0 \end{eqnarray} }$

この方程式の根は $\displaystyle{z = \frac{1}{\beta_1}}$ の一つ。反転可能になるためには $\displaystyle{|z| > 1}$ となればいいので

$\displaystyle{ |z| = \left|\frac{1}{\beta_1}\right| > 1 \\ \\ |\beta_1| < 1 }$

よって $\displaystyle{|\beta_1| \lt 1}$ が条件となります。先程の $\displaystyle{|r| \lt 1}$ の条件とも一致します。

まとめ

MAモデルは同じ自己相関構造をもつ複数のモデルが存在します
反転可能性の制約を入れることで、自己相関構造からMAモデルは一つに定まります
MA(q)の反転可能性条件は以下のMA特性方程式のすべての根が単位円の外側にあることです（ $\displaystyle{|z| > 1}$ ）

$\displaystyle{ \begin{align} \theta(z) &= 1 + \theta_1z + \theta_2z^2 + \dots + \theta_qz^q = 0 \end{align} }$

MAモデルを $\displaystyle{\varepsilon_t}$ を入力とするシステムとして捉える感じが出てきました。制御や信号処理とつながるのは興味深いですね。いつかすべて理解出来る日まで日々精進

2020-08-10

時系列分析：ホワイトノイズとiid

時系列確率過程データ分析 ARMAモデル時系列分析

ホワイトノイズとiid（independent and identically distributed:独立同分布）を混同してしまう事があります。本投稿ではこの2つを整理します。普段の分析ではホワイトノイズとiidをそこまで厳密に区別するケースは少ないと思います。理論を追っていて、前提条件として同一分布を意識するときぐらいでしょうか

ホワイトノイズとiidの違い（まとめ）
- 共通点：
- 違い：
ホワイトノイズの定義
iidの定義
一般的な使い分け
各モデル誤差項の設定
参考文献

ホワイトノイズとiidの違い（まとめ）

共通点：

ホワイトノイズとiidはどちらも確率変数の列がそれぞれ互いに独立です。なので自己共分散は0となります。

違い：

ホワイトノイズはノイズの性質なので平均が0。一方で、iidは独立同分布といっているだけなので、平均値は必ずしも0とは限らない
ホワイトノイズは独立でさえあればいいけど、一方でiidの系列は独立かつ同分布の必要がある

要素を全部まとめると次の表のようになります（意味合い重複した部分あり）

	よく使う表記	平均	自己共分散	独立	同分布	定常性
ホワイトノイズ	$\varepsilon_t \sim \operatorname{W.N.}(\sigma)$	0		必要あり	必要なし	弱定常
iid	$\varepsilon_t \sim \operatorname{iid}(\mu, \sigma)$	$\mu$		必要あり	必要あり	強定常

iidのように同一分布に従うことは強い仮定なので、それに比べて同一分布を仮定しないホワイトノイズは仮定が弱いです。また定常性という視点では、ホワイトノイズは弱定常、iidは強定常になります。

ホワイトノイズの定義

沖本本の定義をそのまま引用するとホワイトノイズの定義は以下のようになっています。ホワイトノイズは平均値0、分散 $\displaystyle{\sigma^ 2}$ 、自己相関（共分散）が0です。

すべての時点 $\displaystyle{t}$ において

$\displaystyle{ \begin{array} \\ E(\varepsilon_t) = 0\\ \gamma_k = E(\varepsilon_t\varepsilon_{t-k}) = \left\{ \begin{array}{ll} \sigma^2,& k = 0 \\ 0, & k\ne 0 \end{array} \right. \end{array} }$

iidの定義

同様にiidの定義は以下です

各時点のデータが互いに独立でかつ同一の分布に従う系列

引用：経済・ファイナンスデータの計量時系列分析 - 沖本竜義 p11

分布が同じで独立な系列であることを要求しているだけなので、そもそもガウシアンのように期待値 $\displaystyle{\mu}$ と分散 $\displaystyle{\sigma^ 2}$ のみで決定される分布かどうかも分かりません。

一般的な使い分け

我々がシミュレーションなどで誤差項を生成する場合、平均値0の正規分布に従うノイズを発生させることが多いと思います。これらのノイズは独立なのでホワイトノイズですが、同一の分布に従うためiidでもあります。いわゆるガウシアンホワイトノイズと言われるもので、iidになります。

各モデル誤差項の設定

沖本本やWikipediaを参考に各モデルの誤差項の仮定を確認してみたいと思います。結論から言うと、MAモデル、ARモデル、ARMAモデルともに撹乱項はホワイトノイズで、iidを仮定していません。ランダムウォークだけなぜかiidを仮定しているみたいです。

	撹乱項の仮定
MAモデル	ホワイトノイズ
ARモデル	ホワイトノイズ
ARMAモデル	ホワイトノイズ
ランダムウォーク	iid

MAモデル

沖本本のMA(q)過程の定義では、MA(q)過程はホワイトノイズの線形和になっています。iidを仮定していません

$\displaystyle{ y_t = \mu + \varepsilon_t + \theta_1\varepsilon_{t-1} + \theta_2\varepsilon_{t-2}+ \dots + \theta_q\varepsilon_{t-q}, \hspace{5mm} \varepsilon_t \sim \operatorname{W.N.}(\sigma^2) }$

引用：経済・ファイナンスデータの計量時系列分析 - 沖本竜義 p24

Wikipediaでも同様に $\displaystyle{\varepsilon_t}$ をwhite noise error termsと記載しています。

ARモデル

ARモデルの撹乱項 $\displaystyle{\varepsilon_t}$ もホワイトノイズです。iidを仮定していません

$\displaystyle{ y_t = c + \phi_1y_{t-1} + \phi_2y_{t-2}+ \dots + \phi_qy_{t-q} + \varepsilon_t, \hspace{5mm} \varepsilon_t \sim \operatorname{W.N.}(\sigma^2) }$

引用：経済・ファイナンスデータの計量時系列分析 - 沖本竜義 p26

Wikipediaでも同様に $\displaystyle{\varepsilon_t}$ をwhite noiseと記載しています。

ARMAモデル

ARMAモデルの撹乱項 $\displaystyle{\varepsilon_t}$ もホワイトノイズです。iidを仮定していません

$\displaystyle{ y_t = c + \phi_1y_{t-1} + \phi_2y_{t-2}+ \dots + \phi_qy_{t-q} + \dots + \varepsilon_t + \theta_1\varepsilon_{t-1} + \theta_2\varepsilon_{t-2}+ \dots + \theta_q\varepsilon_{t-q}, \hspace{5mm} \varepsilon_t \sim \operatorname{W.N.}(\sigma^2) }$

引用：経済・ファイナンスデータの計量時系列分析 - 沖本竜義 p34

ランダムウォーク

ランダムウォークでは撹乱項 $\displaystyle{\varepsilon_t}$ にiidを仮定しています。

$\displaystyle{ y_t = \delta + y_{t-1} + \varepsilon_{t}, \hspace{5mm} \varepsilon_t \sim \operatorname{iid}(0, \sigma^2) }$

引用：経済・ファイナンスデータの計量時系列分析 - 沖本竜義 p106

ランダムウォークをAR(1)過程と考えるならホワイトノイズでも良さそうなものですが、このiidの仮定が重要なのでしょうか。Wikipediaでも独立かつ同分布な Rd 値確率変数族と書いてあります。

なぜランダムウォークだけiidを仮定しているのか分かってません。どなたかご存知でしたら教えて下さい。日々精進

参考文献

経済・ファイナンスデータの計量時系列分析 - 沖本竜義
ランダムウォーク - Wikipedia
Autoregressive–moving-average model - Wikipedia
ホワイトノイズ - Wikipedia
自己回帰移動平均モデル - Wikipedia
独立同分布 - Wikipedia

2020-08-09

時系列分析：後退オペレータについて

データ分析時系列確率過程 ARMAモデル

時系列の勉強をしていると、後退オペレータ（Backshift Operator）というのが出てきます。このオペレータを使って時系列モデルを表現することで、簡潔で見通しの良い議論が可能になるようです。本投稿ではこれについて考察したいと思います。なお後退オペレータは、様々な呼び名があるようで、シフトオペレータ、時間シフトオペレータ、Backshift Operator、Backwardshift Operator など異なる記述で書かれています。

後退オペレータの定義
後退オペレータの多項式表現
後退オペレータの使いどころ
後退オペレータは何なのか
参考文献

後退オペレータの定義

英語版Wikiでは次のように説明してあります

時系列分析において後退オペレータは時系列の要素に作用し、前の要素を生成する

In time series analysis, the lag operator (L) or backshift operator (B) operates on an element of a time series to produce the previous element.

つまり時系列のある時点 $t$ の値 $\displaystyle{X_t}$ を一時点前の $\displaystyle{X_{t-1} }$ に戻すような操作になります。

次のような時系列 $\displaystyle{X}$ があったとすると

$\displaystyle{ X = \{X_1, X_2\, X_3, \dots \} }$

後退オペレータ $\displaystyle{B}$ は以下のように定義されます。

$\displaystyle{ X_{t-1} = BX_t }$

ある時点 $t$ の要素 $\displaystyle{X_t}$ に後退オペレータ $\displaystyle{B}$ を乗ずることで、一時点前の要素 $\displaystyle{X_{t-1}}$ を得ることができます。

同様に二つ前の要素は

$\displaystyle{ \begin{align} X_{t-2} &= BX_{t-1}\\ &= B(BX_{t}) \\ &= B^2X_t \end{align} }$

のようにアクセスできます。一般に $\displaystyle{k}$ 時点前の要素は後退オペレータを用いて次のようにアクセスできます。

$\displaystyle{ X_{t-k} = B^kX_t }$

後退オペレータの多項式表現

後退オペレータを用いて時系列モデル（確率過程）を簡潔に表現するケースは多く見受けられます。その際に次のように後退オペレータ $\displaystyle{B}$ の多項式 $\displaystyle{\phi(B)}$ を使って統一的にモデルを表現します。まずこのことを頭に入れて次の具体例を眺めると良いと思います。

$\displaystyle{ \begin{align} \phi(B) &= \beta_0 + \beta_1B + \beta_2B^2 + \beta_3B^3 + \dots + \beta_kB^k\\ &=\sum_{i=0}^{k}\beta_{i}B^i \end{align} }$

後退オペレータの使いどころ

それでは具体例を見ていきましょう。以下の例はCourseraのここに記載している内容ほぼそのまんまです。

1. ランダムウォーク

まず以下のランダムウォークを考えます。

$\displaystyle{ X_t = X_{t-1} + \varepsilon_t\\ }$

where

$\displaystyle{ \varepsilon_t \sim \operatorname{iid}(0, \sigma^2) }$

ランダムウォークを後退オペレータを用いて表すと

$\displaystyle{ \begin{align} X_t &= X_{t-1} + \varepsilon_t\\ &=BX_{t} + \varepsilon_t\\ \end{align} }$

となります。

さらに後退オペレータの多項式 $\displaystyle{\phi(B)}$ を使って次のようになります。

$\displaystyle{ (1 - B)X_t = \varepsilon_t \\ \phi(B)X_t = \varepsilon_t\\ }$

where

$\displaystyle{ \phi(B) = 1 - B\\ }$

ランダムウォーク自体がシンプルなので、さほど後退オペレータのありがたみは感じません。

2. MA(q)モデル

MA(q)モデルを見てみましょう。ランダムウォークと同じですがMA過程なので $\displaystyle{\varepsilon_t}$ に対して後退オペレータを適用します

$\displaystyle{ \begin{align} X_t &= \varepsilon_t + \beta_1\varepsilon_{t-1} + \beta_2\varepsilon_{t-2} + \dots + \beta_q\varepsilon_{t-q}\\ &= \varepsilon_t + \beta_1B\varepsilon_{t} + \beta_2B^2\varepsilon_{t} + \dots + \beta_qB^q\varepsilon_{t}\\ &= （1 + \beta_1B + \beta_2B^2 + \dots + \beta_qB^q）\varepsilon_{t}\\ \end{align} }$

$\displaystyle{ \begin{align} X_t &=（1 + \beta_1B + \beta_2B^2 + \dots + \beta_qB^q）\varepsilon_{t}\\ &= \phi(B)\varepsilon_{t} \end{align} }$

where

$\displaystyle{ \phi(B) = 1 + \beta_1B + \beta_2B^2\varepsilon_{t} + \dots + \beta_qB^q }$

$\displaystyle{X_t = \phi(B)\varepsilon_t}$ というのは簡潔ですね。

3. AR(p)モデル

またAR(p)過程は、次のようになります

$\displaystyle{ \begin{align} X_t &= \beta_1X_{t-1} + \beta_2X_{t-2} + \dots + \beta_qX_{t-q} + \varepsilon_t\\ &= \beta_1BX_{t} + \beta_2B^2X_{t} + \dots + \beta_qB^qX_{t} + \varepsilon_t\\ \end{align} }$

$\displaystyle{ \begin{align} X_t - \beta_1BX_t + \beta_2B^2X_{t} + \dots + \beta_qB^qX_{t}= \varepsilon_t\\ (1 - \beta_1B + \beta_2B^2 + \dots + \beta_qB^q)X_{t}= \varepsilon_t\\ \phi(B)X_{t}= \varepsilon_t\\ \end{align} }$

where

$\displaystyle{ \phi(B) = 1 - \beta_1B + \beta_2B^2 + \dots + \beta_qB^q }$

$\displaystyle{\phi(B)X_{t}= \varepsilon_t}$ というのも簡潔ですね。

4. ARMAモデル

上記のMA(q)とAR(p)を踏まえてARMAモデルも考えてみます。

$\displaystyle{ X_t = \sum_{i=1}^{p}\alpha_iX_{t-i} + \varepsilon_t + \sum_{i=1}^{q}\beta_i\varepsilon_{t-i}\\ X_t - \sum_{i=1}^{p}\alpha_iB^iX_{t} = \varepsilon_t + \sum_{i=1}^{q}\beta_iB^i\varepsilon_{t}\\ (1 - \sum_{i=1}^{p}\alpha_iB^i)X_{t} = (1 + \sum_{i=1}^{q}\beta_iB^i)\varepsilon_{t}\\ \phi_{\alpha}(B)X_t = \phi_{\beta}(B)\varepsilon_t\\ }$

where

$\displaystyle{ \phi_{\alpha}(B) = 1 - \sum_{i=1}^{p}\alpha_iB^i\\ \phi_{\beta}(B) = 1 + \sum_{i=1}^{q}\beta_iB^i\\ }$

ARMAモデルは $\displaystyle{\phi_{\alpha}(B)X_t = \phi_{\beta}(B)\varepsilon_t}$ のように大変スッキリした形で記述できることがわかります。ARモデルとMAモデルの表現も統一感があります。

5. 差分作用素

差分作用素は後退オペレータでも記述できます

１階差の場合:

$\displaystyle{ \begin{align} \nabla X_t &= X_t - X_{t-1}\\ &= (1 - B)X_t \end{align} }$

2階差の場合:

$\displaystyle{ \begin{align} \nabla^2 X_t &=(X_t - X_{t-1}) - (X_{t-1} - X_{t-2}) \\ &= (1 - B)X_t - (1 - B)X_{t-1}\\ &= (1 - B)(X_t - X_{t-1})\\ &= (1 - B)^2X_t\\ \end{align} }$

一般に

$\displaystyle{ \nabla^iX_t = (1 - B)^iX_t }$

後退オペレータは何なのか

ところで後退オペレータの実態は何なのでしょうか。ある関数 $\displaystyle{f(t)}$ に乗ずることで、 $\displaystyle{t}$ をシフトさせるような作用素といえば、複素数が思い浮かびます。厳密ではないですが、後退オペレータは複素数そのものと考えて問題なさそうな気がします。

たとえば関数 $\displaystyle{ x(t)}$ があるとして、フーリエ級数に展開できるとすると

$\displaystyle{ x(t) ＝ \sum_{k=\infty}^{\infty}c_k\mathrm{e}^{jkt} }$

等間隔 $\displaystyle{d}$ でサンプリングされたすると

$\displaystyle{ X(n) = \sum_{k=\infty}^{\infty}c_n\mathrm{e}^{jkdn} }$

これを $\displaystyle{X(n-1)}$ にしたい場合、 $\displaystyle{\mathrm{e}^ {-jdn}}$ を乗じてあげれば、なんかよさそうです。

$\displaystyle{ \begin{align} \mathrm{e}^{-jdn}X(n) &= \mathrm{e}^{-jdn}\sum_{k=\infty}^{\infty}c_n\mathrm{e}^{jkdn}\\ &= \sum_{k=\infty}^{\infty}c_n\mathrm{e}^{j(k-1)dn}\\ &= X(n - 1) \end{align} }$

なので、この場合 $\displaystyle{B = \mathrm{e}^ {-jdn}}$ で問題ないかなと思うのですが、どなたかご存知でしたら教えて下さい。

シフト作用素 - Wikipedia

追記：関数解析学の分野でシフト作用素として一般化されているものでした。折を見て、整理し直したいと思います。日々精進

参考文献

Lag operator | WIKIPEDIA
経済・データの計量時系列分析 (統計ライブラリー) – 沖本竜義
Practical Time Series Analysis | Coursera
Appendix B COMPLEX VARIABLES, THE SPECTRUM, AND LAG OPERATOR
時系列解析入門 – 北川源四郎

2020-08-02

はてなブログで日々数式を含んだ記事を書く

Python TeX ブログ正規表現

はてなブログを復活して2週間。はてなブログでTeXを使って数式を書くのが大変で当初ずっとイライラしていた。現時点での私の対応策を書きたいと思う。

まずQiitaと比べて圧倒的に面倒な点は次：

下書きがリアルタイムでレンダリングできない（数式、リンクを含む一部）
TeXの記述が独特で特にエスケープがめちゃくちゃ面倒

1の対策. Typoraで下書き

[1]に関しては、Typora — a markdown editor, markdown reader.を使うことで、markdown + TeX をリアルタイムにて下書きを確認している。私はTeXもそんなに慣れているわけではないのでリアルタイムにレンダリングしてくれないと苦しい

Typoraで下書きをすれば、TeXもリアルタイムに確認できる

f:id:messefor:20200802190337p:plain — ypora

TeXはmathpixで

ちなみにTeXコードも複雑な式はMathpix Snipを使って、手書きノートの写真から直接TeXコードに変換している。

f:id:messefor:20200802191311p:plain — mathpix

たとえばこれが、

f:id:messefor:20200802194010p:plain — formula

mathpix で読み込むとこうなる。 f:id:messefor:20200802194307p:plain

もちろん完璧ではないが、精度も十分。めちゃ便利。なおこれはpngで出力したもの。

2の対策. スクリプトで置換

当初、はてなブログでのTeX エスケープの方法を探り探りやっていたが、結局以下の記事に記載してあることが全てだった。

ano3.hatenablog.com

上記ブログに設置してあるフォームから変換できるのだが、作業の流れ上Typoraで書いた.mdファイルをそのまま一気に変換したかったので、とりいそぎ自前でPythonでスクリプトを書いた。

実行部分

path_art = 'article.md'
path_art_cnvd = 'article_new.md'

# 変換対象のテキストファイルの読み込み
with open(path_art, 'r') as f:
    artile_strings = f.read()

# 変換 
rht = ConvertHatenaTex()
artile_strings_new = rht.replace(artile_strings)

# 変換後のファイルを保存
with open('article_new.md', 'w') as f:
    f.write(artile_strings_new)

やっていることは、$$~$$と $~$ で囲まれたブロックを見つけてきて、愚直に必要部分を置換しているだけ。

全体

"""Convertace article with TEX syntax with Hatena Blog Syntax """
import os
import re
from typing import List, Callable
from copy import deepcopy


def replace_m(x: str, rep_map: List[dict]) -> str:
    for pattern, repl in rep_map:
        x = x.replace(pattern, repl)
    return x

def add_tags(x: str, tag_pre: str, tag_post: str) -> str:
    """Concatenate string and tags"""
    return tag_pre + x + tag_post


class ConvertHatenaTex:
    """Convert TEX syntax to Hatena Blog TEX Syntax.

    """

    def __init__(self):

        self.pattern_block = r'\$\$([\S\s]*?)\$\$'
        self.pattern_inline = r'\$([a-zA-Z0-9!-/:-@¥[-`{-~ ]+)\$'

        self.rep_map_block = [('[', r'\['), (']', r'\]'),]
        self.rep_map_inline = [('[', r'\\['), (']', r'\\]'),
                                ('_', r'\_'), ('^', '^ '),
                                (r'\{', r'\\{') ,(r'\}', r'\\}')]


        tag_pre_b = '<div  align="center">\n[tex:\displaystyle{\n\n'
        tag_post_b = '\n}]\n</div>\n'

        # tag_pre_b = '<pre>\n[tex:\n\n'
        # tag_post_b = '\n]\n</pre>\n'

        self.tags_block = (tag_pre_b, tag_post_b)

        tag_pre_i = '[tex:\displaystyle{'
        tag_post_i = '}]'

        # tag_pre_i = '[tex:'
        # tag_post_i = ']'

        self.tags_inline = (tag_pre_i, tag_post_i)

    def replace(self, artile_strings: str) -> str:
        x = deepcopy(artile_strings)
        x = self._replace(x, self.pattern_block, self.proc_in_blocks)
        x = self._replace(x, self.pattern_inline, self.proc_in_lines)
        return x

    def proc_in_blocks(self, x: str) -> str:
        """Process replacing lines in blocks"""
        x = replace_m(x, self.rep_map_block)
        x = add_tags(x, *self.tags_block)
        return x

    def proc_in_lines(self, x: str) -> str:
        """Process replacing in lines"""
        x = replace_m(x, self.rep_map_inline)
        x = add_tags(x, *self.tags_inline)
        return x

    def _replace(self, artile_strings: str, pattern: str,
                            rep_func: Callable[[str], str]) -> str:
        """Extract block/inline parts and replace"""
        prog = re.compile(pattern)
        results = prog.finditer(artile_strings)

        replace_strs = []
        n_init = 0
        for mg in results:

            str_chunk = mg.group(1)

            str_chunk = rep_func(str_chunk)

            n_start, n_end = mg.span()
            str_pre = artile_strings[n_init:n_start]

            replace_strs.append(str_pre)
            replace_strs.append(str_chunk)
            n_init = n_end

        replace_strs.append(artile_strings[n_end:])
        artile_strings_new = ''.join(replace_strs)

        return artile_strings_new



if __name__ == '__main__':

    path_art = 'article.md'
    path_art_cnvd = 'article_new.md'

    # Load a article to be replaced
    with open(path_art, 'r') as f:
        artile_strings = f.read()

    rht = ConvertHatenaTex()
    artile_strings_new = rht.replace(artile_strings)

    # Save the replaced article
    with open(path_art_cnvd, 'w') as f:
        f.write(artile_strings_new)

人生残り短い中年に残された時間は少ない。日々のアウトプットにあまり時間を使いすぎるのは精神衛生上良くない。

オススメの方法があったら教えてください。