フリーの数学ソフトMAXIMAの力を借りて差分方程式の安定条件、ARMA(1,1)モデルの反転可能性、因果性等を少し探ってみた。末尾に示したようなEndersのApplied Econometric Time seriesの邦訳書が出ているので、そこで取り上げられているARMA(1.1)の計算例を中心にして、必要に応じてブロックウェル/デービス(Brockwell and Davis)のテキストも参照しながら探ってみる。
・
最初に差分方程式のパッケージsolve_rec を読み込んでおく。


安定条件について幾つかの差分方程式を解きながら検討してみる。

y(t)=0.2y(t-1)+0.35y(t-2)
のケース

この差分方程式を解いてみる。



%k1　%k2　は任意定数をあらわす。k1=1 k2=1 として特殊解を求めると下記のような式を得る。



この式を使ってグラフを描く。時間とともに急減している。



y(t)=0.7y(t-1)+0.35y(t-2)のケース

差分方程式を解いてみる



%k1　%k2　は任意定数をあらわす。k1=1 k2=1
として特殊解を求めると下記のような式を得る



この式を使ってグラフを描く。時間とともに発散している。



y(t)=1.6y(t-1)-0.9y(t-2)
のケース
2階差分方程式を解いてみる



この式を使ってグラフを描く。時間ととも振幅しながら収束している。

特性方程式と反転特性方程式

characteristic equation and inverse characterisitic equation

安定性を確かめるため特性方程式(characteristic equation)を解いてみる。



特性根(characteristic roots)の絶対値(absolute value)はいづれも１よりも小さい。単位円の内側(with in the unit circle)にあり安定条件を満たしていることが分かる。

この複素数の偏角を計算し三角関数を使ってグラフを描いてみる。



上記の差分方程式についてLag演算子を使ってラグ多項式を作り、その多項式から下記のような反転特性方程式(inverse characteristic equation)を立てる。
　1-1.6z+0.9z²=0　
その根を求めてみる。



特性根zの絶対値はすべて１よりも大きい。つまり単位円の外側(without the unit circle)にある.反転特性方程式を使った場合には単位円の外側にあることが安定条件となる。　
２つの特性方程式の根の絶対値0.9486と1.0540は互いに逆数の関係にある。(these characteristic roots are reciprocals each other)
時系列分析では特性方程式の立て方で安定条件は単位円の内側あるいは外側にあることが必要となり説明が変わるので話が紛らわしい。
　ささやかな経験ではラグ多項式を使ってＡＲＭＡモデルを説明するテキストが多いような気もするのでラグ多項式から反転特性方程式を立てて、すべての特性根の絶対値が１より大きい、つまり単位円の外側に存在することが安定条件を満たすと理解しておく方が便利かも知れない。

 2階差分方程式の根と係数の関係
stability conditions of 2nd order difference equation

y(t)-a*y(t-1)-b*y(t-2)＝０の形となる２階差分方程式の係数と安定条件の関係（根と係数の関係）を探ってみる。



２階差分方程式の係数a,bが上記のグラフの三角形部分に存在していれば制約条件を満たし、つまり安定条件が満たされる。

y(t)-0.7y(t-1)-0.35y(t-2)=0　の差分方程式で検討するとa=0.7 b=0.35 は何れも１以下であるが　b=1-a=1-0.7=0.3 から　bは0.3より小さいことが必要だが　b=0.35 で0.3を超えてしまい三角形の外にあり必要な条件を満たしていない。
念のためラグ多項式に基づく反転特性方程式を
解いてみると(1-0.7L-0.35L²)y(t)から
1-0.7z-0.35z²=0 を解くと


全ての根が単位円の外側にある要件を満たしておらず１つの根は単位円内にあるので安定条件を満たさずグラフが示すように発散している。

特性方程式を
z²-0.7z-0.35=0 の形で解くと
全ての根が単位円の内側にあることが安定条件となるが、それを満たしているか調べてみる。



全ての根が単位円の内側にある要件を満たしていない。
このように特性方程式の立て方次第で根は異なってくるが、相互に逆数の関係にある。従って単位円の内側とか外側とか説明が異なってくる。

ARMA(1,1)モデルの検討

因果性(causality),反転可能性(invertibility)の検討

ARMA(1,1)過程をラグ演算子を使ってAR項とMA項を整理してみる。ここでa,bは係数,e ホワイトノイズw₁,w₂...w_j は係数　とする。
y(t)=a*y(t-1)+e(t)+b*e(t-1)
(1-aL)y(t)=(1+bL)e(t)　
ここで多項式(1-aL)、(1+bL)が共通の根を持たず、それぞれ全ての根が単位円の外側にある場合、つまり|z|>1 であることが、因果性、反転可能性を満たすとともに多項式の分数
(1+bL)/(1-aL)　がw₀+w₁L+w₂L²+w₃L³+といった無限の項からなる多項式で表せることが証明されている。(数学的に正確な説明は末尾の参考文献などを參照されたい。) そこで未定係数法を使ってw₀,w₁,w₂...を計算することによりMA(∞),あるいはAR(∞)の形で表現ができることが確認できる。以下ではARMA(1,1)が必要条件を満たせば、MA(∞)で表現できることを未定係数法によって確かめている。






上記と同様の方法で未定係数法を使い、(1-aL)/(1+bL)の新しい多項式の係数を求めて、ARMA(1,1)がAR(∞)で表現できることも確かめられる。
 
次のモデルはARMA(1,1)の形をしている。
y(t)=0.7y(t-1)+e(t)-0.7*e(t-1) これをラグ多項式で整理する
(1-0.7L)y(t)=(1-0.7L)e(t) 
結局 y(t)がホワイトノイズであることつまり y(t)=e(t) を回りくど表現したに過ぎないことが分かる .
しかもAR項とMA項の反転特性方程式を解くと、それぞれ
1-0.7z=0 から共通の根10/7≒1.4285が存在することになりARMA過程を定義できない。
そのため、多項式(1-aL)、(1+bL)が共通の根を持たず、それぞれ全ての根が単位円の外側にあることがARMA過程の存在する条件となっている。

因果性(causality)について

AR(1)過程でy(t)について反復代入の計算をすると



係数aが |a|<1 であればy(t)はMA(∞)で表現できる。つまり現在の値y(t)は現在及び過去のホワイトノイズ（ショック）に依存していることを表している。これは現実感として納得しやすい。

もし|a|>1の時には数式を少し操作して
y(t+1)=a*y(t)+e(t+1) で表し、両辺をaで除して下記のように整理すると定常解が得られる。




定常解を得るがy(t)は未来のホワイトノイズeに依存することになり 不自然な解となる。そこで現実感のあるARモデル、つまり現在及び過去の 値にだけ依存する時系列モデルを因果的(causal)として区別するため 因果性(causality)の条件を課すことが行われている。その条件は AR(p)過程がラグ多項式a(L)=1−a₁L-a₂L^2-∙ ∙∙--a_pL^p  において Lをｚに置き換えて得る反転特性方程式a(z)=1−a₁z-a₂z^2-∙ ∙∙--a_pz^p の全ての根が単位円の外側にあること、つまり|z|>1  であることが 必要とされる。 これと同様にMA(q)過程でMA項のラグ多項式b(L)=1−b₁L-b₂L²-∙ ∙∙--b_pL^p
b(z)=1−b₁z-b₂z²-∙ ∙∙--b_pz^p の全ての根が単位円の外側にあること、つまり|z|>1  であるときには反転可能(invertible)と言われる。つまりMA項がAR(∞) で表現できる。正確な数学的な説明は末尾の参考文献を参照されたい。


次のようなARMA(1,1)モデルを検討してみる。

y(t)=-0.7y(t-1)+e(t)-0.7e(t-1)

これをラグ演算子を使って式を整理すると

(1+0.7L)y(t)=(1-0.7L)e(t)

AR項の反転特性方程式を解くと


zの絶対値は１より大きいので安定的、因果的である。

MA項の反転特性方程式を解くと



zの絶対値は１より大きいので反転可能(invertible)である。

このARMA(1,1)モデルは因果的であり、かつ、反転可能なモデルだと分かる。

ARMA(1,1)の自己相関関数(ACF)

自己相関関数を検討するため、以下のような表記を使うことにする。



ARMA(1,1)の自己相関の計算過程は煩雑になるが以下のようになる。











γ(0)とγ(1)が求まればk=1の自己相関ρ(k=1)はγ(1)/γ(0)＝ρ(1)で計算できる。以下では具体的な計算で確かめてみる。


ARMA(1,1)の自己相関関数 rho を以下のように定義すると rho(0)=1


t=0から8までの自己相関係数を求める。



グラフにすれば




自己相関関数は
ρ(k)=a*ρ(k-1)といった1階差分方程式で表せるので初期条件を与えて解くと



この式を利用して関数を作ると


このように同様の結果を得る。

数多くの時系列分析のテキストが出版されているが、翻訳書の場合は訳注とか補足的説明が加えられていることがあり、原文よりも分かりやすいといったメリットも感じられる。ブッロックウエル／デービスのテキストは入門というタイトルがついているが、個人的にはかなり難しく手ごわいという感じがする。Rが使用できる環境であればwww.r-project.orgからpackage  itsmr をインストールすることでITSM2000と同様のソフトが使えるので学習には便利だと思う。


参考文献
ウォルター・エンダース著 新谷元嗣、藪友良 訳
(2019)実証のための計量時系列分析 有斐閣(Walter Enders Applied Econometric time series Wiley & Sons Inc)


ブロックウェル/ デービス 著
逸見功,田中稔,宇佐美嘉弘,渡辺則生 訳
(2004)入門時系列解析と予測（改訂第２版）
シーエーピー出版株式会社(Brockwell and Davis Introduction to time series and forecasting Springer-Verlag)