statsmodelsの季節分解で実装されているアルゴリズム

前回に続いて、statsmodelsの季節分解の話です。
参考(前回の記事): statsmodelsを利用した時系列データの季節分解のやり方

前回の記事は使い方でしたが、今回はソースコードを参照しながらどのような計算方法で季節分解が実装されているのかを見ていきます。
ちなみに、僕の環境のバージョンは以下の通りです。将来のバージョンでは仕様が変わる可能性もあるのでご注意ください。

$ pip freeze | grep statsmodels
> statsmodels==0.13.2

ドキュメントにソースコードが載ってるページがあるので、そこを参照します。
ソース: statsmodels.tsa.seasonal — statsmodels

ソースの先頭に以下のようにコメントが書かれている通り、移動平均を使って実装されています。

"""
Seasonal Decomposition by Moving Averages
"""

利用するデータですが、前回の記事と同じ、このブログのPVです。

# データ件数
print(len(df))
# 140

# 2週間分のデータ
print(df.tail(14))
"""
              pv
date
2022-10-24  2022
2022-10-25  2140
2022-10-26  2150
2022-10-27  1983
2022-10-28  1783
2022-10-29   847
2022-10-30   793
2022-10-31  1991
2022-11-01  2104
2022-11-02  1939
2022-11-03  1022
2022-11-04  1788
2022-11-05   830
2022-11-06   910
"""

それでは順番にやっていきましょう。モデルは加法モデルの方を取り上げます。乗法モデルもトレンド成分を抽出するところまでは一致していますし、その後の処理から異なりますが加法と乗法の違いを考えれば普通に理解できると思います。

そして、ここが重要なのですが、まず周期が奇数の場合を見ていきます。今回は1週間である7日です。なんでこれが重要かというと、偶数の場合は少し挙動が特殊だからです。

季節分解では、元の時系列データからトレンド成分と季節成分を取り出します(残りが残差)が、処理の順番もこの順番になっていて、最初にトレンド成分、次に季節成分を取り出すようになっています。これ順番逆にすると結果変わりますし、それぞれメリットデメリットあるのですがトレンド成分を先にFIXすることを選ばれているようです。

では、トレンド成分の抽出を見ていきましょう。これ、”by Moving Averages”のコメント通り、移動平均です。two_sided=True (デフォルト)の設定場の場合、その対象日を挟むようにして、前後の期間から取った移動平均をその日のトレンド成分の値とします。

7日だったら、3日前,2日前,1日前,当日,1日後,2日後,3日後 の平均を使います。

pandasで計算するなら、rolling()で移動平均取って、shift()でずらすと同じ値が得られます。
モデルで取り出したトレンド成分と、pandasで自分で計算した値見てみましょう。

# モデルで計算
from statsmodels import api as sm


decompose_result = sm.tsa.seasonal_decompose(
    df["pv"],
    period=7,
)
print(decompose_result.trend)
"""
date
2022-06-20            NaN
2022-06-21            NaN
2022-06-22            NaN
2022-06-23    1409.857143
2022-06-24    1408.714286
                 ...     
2022-11-02    1495.285714
2022-11-03    1512.000000
2022-11-04            NaN
2022-11-05            NaN
2022-11-06            NaN
Name: trend, Length: 140, dtype: float64
"""

# pandasで計算。7日移動平均をとって、3日シフトする
print(df["pv"].rolling(7).mean().shift(-3))
"""
date
2022-06-20            NaN
2022-06-21            NaN
2022-06-22            NaN
2022-06-23    1409.857143
2022-06-24    1408.714286
                 ...     
2022-11-02    1495.285714
2022-11-03    1512.000000
2022-11-04            NaN
2022-11-05            NaN
2022-11-06            NaN
Name: pv, Length: 140, dtype: float64
"""

データの大部分…て略されてますが、これらは一致します。(極小値の誤差はあり得ますが)

two_sided = False とすると、前後のデータではなくその日含めた過去のデータだけ使われるようになるので、6日前〜当日のデータでの移動平均になります。

比較できるように、最初の2個の値が見れるようにprintしました。先頭のNaNの数が6個になってますね。表示てませんが、代わりに末尾のNaNは無くなります。

decompose_result = sm.tsa.seasonal_decompose(
    df["pv"],
    period=7,
    two_sided=False
)
print(decompose_result.trend.iloc[: 8])
"""
date
2022-06-20            NaN
2022-06-21            NaN
2022-06-22            NaN
2022-06-23            NaN
2022-06-24            NaN
2022-06-25            NaN
2022-06-26    1409.857143
2022-06-27    1408.714286
Name: trend, dtype: float64
"""

さて、トレンド成分の計算方法はわかったので、季節成分の話に戻ります。モデルも two_sided=Trueの最初の例の方に話戻してこちらで進めます。(この記事真似して再現する人は上のtwo_sided=Falseのコードブロックをスキップして実行してください)

季節成分は、トレンド成分を取り除き終わったデータから計算します。(最初と最後の数件のデータはNaNになっているのでこれらは使いません。)
説明が難しいのですが、今回みたいに曜日ごとの周期であれば、月曜の平均、火曜の平均、水曜の平均と、周期分の平均値を取り出し、さらにこうして計算した平均値たちからその平均を引いて標準化します。

言葉で書くとわかりにくいのでコードでやってみます。

import numpy as np


# 先述のトレンド成分の計算
df["pv_trend"] = df["pv"].rolling(7).mean().shift(-3)

# トレンド成分を取り除く
df["pv_detrended"] = df["pv"] - df["pv_trend"]

# 曜日ごとの平均を取る
period_averages = np.array([df["pv_detrended"][i::7].mean() for i in range(7)])
print(period_averages)
"""
[ 243.2556391   397.69172932  341.57894737  257.54285714  139.98496241
 -720.2481203  -675.17293233]
"""

# 曜日ごとの平均から、さらにそれらの平均を引く
period_averages -= period_averages.mean()
print(period_averages)
"""
[ 245.450913    399.88700322  343.77422127  259.73813104  142.18023631
 -718.0528464  -672.97765843]
"""

# 比較用。statsmodelsが算出した季節成分
decompose_result = sm.tsa.seasonal_decompose(
    df["pv"],
    period=7,
)
print(decompose_result.seasonal[: 7])
"""
date
2022-06-20    245.450913
2022-06-21    399.887003
2022-06-22    343.774221
2022-06-23    259.738131
2022-06-24    142.180236
2022-06-25   -718.052846
2022-06-26   -672.977658
Name: seasonal, dtype: float64
"""

データの形式が違いますが、モデルの結果と値が一致してるのがわかりますね。

以上で、トレンド成分と季節成分が計算できました。あと、statsmodelsは残差を計算できますがこれは単にトレンド成分と周期成分を元のデータから取り除いてるだけです。 df[“pv_detrended”] はトレンド成分除去済みなのでここから季節成分も引いてみます。

# np.tile で同じ値を繰り返す配列を作って引く
print(df["pv_detrended"] - np.tile(period_averages, 20))
"""
date
2022-06-20           NaN
2022-06-21           NaN
2022-06-22           NaN
2022-06-23     74.404726
2022-06-24     36.105478
                 ...    
2022-11-02     99.940064
2022-11-03   -749.738131
2022-11-04           NaN
2022-11-05           NaN
2022-11-06           NaN
Name: pv_detrended, Length: 140, dtype: float64
"""

# モデルの残差
print(decompose_result.resid)
"""
date
2022-06-20           NaN
2022-06-21           NaN
2022-06-22           NaN
2022-06-23     74.404726
2022-06-24     36.105478
                 ...    
2022-11-02     99.940064
2022-11-03   -749.738131
2022-11-04           NaN
2022-11-05           NaN
2022-11-06           NaN
Name: resid, Length: 140, dtype: float64
"""

以上が、statsmodelsにおける seasonal_decompose の実装の説明になります。

さて、ここからはちょっと補足で、周期が偶数の場合の挙動になります。月次データだと12を使うことが多いので周期が偶数というのはよくあるケースです。
説明をコンパクトにするために、周期(period=)4を例に取り上げます。

周期7の時、3日前から3日後までのデータの平均でトレンド成分を取り出してましたが、こういう期間の取り方ができるのは周期が奇数だったからで、偶数だとこうはいきません。

では、どうするのかというと、例えば周期が4日だったら、前日から翌日までの3日分のデータはそのまま使い、2日前と2日後のデータをそれぞれ0.5倍したものを計算に入れます。

# 周期4の場合のトレンド成分
decompose_result = sm.tsa.seasonal_decompose(
    df["pv"],
    period=4,
)
print(decompose_result.trend.iloc[: 5])
"""
date
2022-06-20         NaN
2022-06-21         NaN
2022-06-22    1719.750
2022-06-23    1567.250
2022-06-24    1299.125
Name: trend, dtype: float64
"""

# 以下、 2022-06-22    1719.750 が算出されるまでの計算式
# 元のデータ
df["pv"].iloc[: 5]
"""
date
2022-06-20    1703
2022-06-21    1758
2022-06-22    1732
2022-06-23    1744
2022-06-24    1587
Name: pv, dtype: int64
"""

# 計算 2日前〜2日後の5日分のデータから算出。ただし、最初と最後の日はウェイトが0.5
(1703*0.5 + 1758 + 1732 + 1744 + 1587*0.5)/4
# 1719.75

rolling()とshift()では同じ値を得られなかったので戸惑ったこともありましたが、まぁ、仕掛けがわかってしまえば簡単ですね。その日を挟んだ4日分のデータを考慮する手段としても一定の妥当性があると思います。

で、ここからが僕的には納得がいってないのですが、two_sided=Falseで周期を偶数(今の例では4)にした場合の挙動は少し不思議です。その日を含む過去4日分のデータを使えるので、単純に4日間の平均を取ればいいのに、そういう実装になっておらず、two_sided=Trueの場合の結果を単純にshiftしたものを使っています。要するに4日前の0.5倍と3日前から1日前、そして当日の0.5倍のデータを使ってます。two_sided=Falseにしたコードが以下ですが、上のTrue(省略した場合の値)と値が一致しているのがわかると思います。

decompose_result = sm.tsa.seasonal_decompose(
    df["pv"],
    period=4,
    two_sided=False,
)
print(decompose_result.trend.iloc[: 7])
"""
date
2022-06-20         NaN
2022-06-21         NaN
2022-06-22         NaN
2022-06-23         NaN
2022-06-24    1719.750
2022-06-25    1567.250
2022-06-26    1299.125
Name: trend, dtype: float64
"""

ここの挙動だけは将来のバージョンで修正されるのではと思っているのですが、
一旦今はこういう作りになっているということを頭に置いた上で気をつけて使うしかないようです。

statsmodelsを利用した時系列データの季節分解のやり方

要するに、seasonal_decomposeメソッドの使い方の紹介記事です。これもとっくに書いたと思っていたら書いてなかったのでまとめておきます。この記事では季節分解の概要の説明とライブラリの使い方を紹介します。そして、これの次の記事でstatsumodelsがどのような実装で季節分解を行っているのかを解説する予定です。
参考: statsmodels.tsa.seasonal.seasonal_decompose — statsmodels

現実の時系列データは、何かしらの季節性を持っていることが多くあります。季節って単語で言うと、春夏秋冬や、何月、といった粒度のものが想像されやすいですが、1週間の中で見ても曜日の傾向とか、1日の中でも時間帯別の違いなどがあります。

時系列データからこの季節に依存する部分を取り出し、季節成分と、トレンド成分、そして残差へと分解する手法が今回紹介する季節分解です。Wikipediaでは季節調整と書いてあります。また、基本成分など、別の用語を使ってる人もいるようです。(どれが一番メジャーなんだろう。)

定式化しておくと、元の時系列データ$Y_{t}$をトレンド成分$T_t$と季節成分$S_t$、そして残差$e_t$を用いて、
$$
Y_t = T_t + S_t + e_t
$$
と分解することを目指します。上記のは加法モデルと呼ばれる形で、和の代わりに積で分解する乗法モデル、
$$
Y_t = T_t * S_t * e_t
$$
もあります。

季節成分$S_t$は周期性を持っているので、その周期を$p$とすると、$S_{t}=S_{t+p}$を満たします。

具体的に例を見るのが早いと思うので、やっていきましょう。サンプルとして用意したデータはこのブログのpv数です。インデックスを日付けにしていますが、こうしておくとライブラリでplotしたときにx軸に日付が表示されるで便利です。ただ、通常の通し番号のindexでも動きます。

# データ件数
print(len(df))
# 140

# 2週間分のデータ
print(df.tail(14))
"""
              pv
date            
2022-10-24  2022
2022-10-25  2140
2022-10-26  2150
2022-10-27  1983
2022-10-28  1783
2022-10-29   847
2022-10-30   793
2022-10-31  1991
2022-11-01  2104
2022-11-02  1939
2022-11-03  1022
2022-11-04  1788
2022-11-05   830
2022-11-06   910
"""

最後に元のデータのグラフも出るので可視化しませんが、上の例見ても10/29,30や11/5,6など土日にpvが減っていて逆に平日多く、曜日ごと、つまり7日周期がありそうな想像がつきます。詳細省略しますが、自己相関等で評価してもはっきりとその傾向が出ます。

では、やってみましょう。まず、分解自体はライブラリにデータを渡して周期を指定するだけです。

from statsmodels import api as sm


decompose_result = sm.tsa.seasonal_decompose(
    df["pv"],
    period=7,  # 周期を指定する
)

結果は以下のプロパティに格納されています。DecomposeResultというデータ型で、ドキュメントはこちらです。
参考: statsmodels.tsa.seasonal.DecomposeResult — statsmodels

順番に表示していきます。

# データの数。
decompose_result.nobs[0]
# 140

# 元のデータ
print(decompose_result.observed[: 10])
"""
date
2022-06-20    1703.0
2022-06-21    1758.0
2022-06-22    1732.0
2022-06-23    1744.0
2022-06-24    1587.0
2022-06-25     654.0
2022-06-26     691.0
2022-06-27    1695.0
2022-06-28    1740.0
2022-06-29    1655.0
Name: pv, dtype: float64
"""
# トレンド成分
print(decompose_result.trend[: 10])
"""
date
2022-06-20            NaN
2022-06-21            NaN
2022-06-22            NaN
2022-06-23    1409.857143
2022-06-24    1408.714286
2022-06-25    1406.142857
2022-06-26    1395.142857
2022-06-27    1370.714286
2022-06-28    1345.285714
2022-06-29    1347.142857
Name: trend, dtype: float64
"""

# 季節成分
print(decompose_result.seasonal[: 14])
"""
date
2022-06-20    245.450913
2022-06-21    399.887003
2022-06-22    343.774221
2022-06-23    259.738131
2022-06-24    142.180236
2022-06-25   -718.052846
2022-06-26   -672.977658
2022-06-27    245.450913
2022-06-28    399.887003
2022-06-29    343.774221
2022-06-30    259.738131
2022-07-01    142.180236
2022-07-02   -718.052846
2022-07-03   -672.977658
Name: seasonal, dtype: float64
"""
# 残差
print(decompose_result.resid[: 10])
"""
date
2022-06-20          NaN
2022-06-21          NaN
2022-06-22          NaN
2022-06-23    74.404726
2022-06-24    36.105478
2022-06-25   -34.090011
2022-06-26   -31.165199
2022-06-27    78.834801
2022-06-28    -5.172718
2022-06-29   -35.917078
Name: resid, dtype: float64
"""

トレンド成分が最初の3項NaNになっているのは、アルゴリズムの都合によるものです。その日を中心とする前後で合計7日(周期分)のデータで移動平均をとっており、要するに、過去の3日、当日、次の3日間、の合計7日分のデータがないと計算できないので最初の3日と、表示していませんが最後の3日間はNaNになっています。この辺の挙動は推定時の引数で調整できます。

次に季節成分は14日分printしましたが、最初の7日間の値が繰り返されて次の7日間でも表示されているのがわかると思います。ずっとこの繰り返しです。

残差は元のデータからトレンド成分と季節成分を引いたものになります。トレンド成分や季節成分に比べて値が小さくなっていて、今回のデータではトレンドと季節である程度分解が綺麗に行えたと考えられます。

さて、データが取れたのでこれを使えばmatplotlib等で可視化できるのですが、大変ありがたいことにこのDecomposeResultが可視化のメソッドを持っています。
少し不便なところは、そのplotメソッドがfigsizeとかaxとかの引数を受け取ってくれないので、微調整とかしにくいのですよね。
個人的な感想ですが、デフォルトでは少しグラフが小さいのでrcParamsを事前に変更してデフォルトのfigsizeを大きくし、それで可視化します。

import matplotlib.pyplot as plt
#  注: DecomposeResult.plot()を実行するだけなら matplotlibのimportは不要。
#       今回画像サイズを変えるためにインポートする

# figure.figsizeの元の設定を退避しておく
figsize_backup = plt.rcParams["figure.figsize"]

# 少し大きめの値を設定
plt.rcParams["figure.figsize"] = [10, 8] 

# 可視化
decompose_result.plot()
plt.show()

# 設定を元に戻す
plt.rcParams["figure.figsize"] = figsize_backup

これで出力される画像が以下です。

一番上が元のpvです。冒頭に書いた通り、生データそのまま見ても周期性が明らかですね。

トレンド成分を見ていくと、お盆の時期にアクセスが減っていますが、その後順調にアクセスが伸びていることがわかりますね。

残差で見ると大きくアクセスが減っているのはそれぞれ祝日に対応しています。
8月に1日だけ異常にアクセス伸びた日がありますがこれは謎です。

説明や例をいろいろ書いてきたので長くなりましたが、基本的には、seasonal_decompose()で分解して、plot()で可視化してそれで完成という超お手軽ライブラリなので、時系列データが手元にあったらとりあえず試す、くらいの温度感で使っていけると思います。

最後に、seasonal_decompose にオプション的な引数が複数あるので使いそうなものを説明してきます。

まず、model= は、 “additive”,”multiplicative” の一方をとり、加法的なモデルか乗法的なモデルを切り替えることができます。デフォルトは、”additive”です。

two_sided= は True/Falseの値をとり、デフォルトはTrueです。これはトレンド成分の抽出方法を指定するもので、Trueの時は、その日を挟むように前後の日付から抽出しますが、Falseの場合は、その日以前の値から算出されます。True or False で並行移動するイメージです。
次回の記事で詳細書こうと思いますが、周期が偶数か奇数かで微妙に異なる挙動をするので注意が必要です。

extrapolate_trend= はトレンド成分の最初と最後の欠損値を補完するための引数です。1以上の値を渡しておくと、その件数のデータを使って最小二乗法を使って線形回帰してトレンドを延長し、NaNをなくしてくれます。使う場合はその回帰が妥当かどうか慎重にみて使う必要がありそうです。

Python(statsmodels)でインパルス応答関数の計算

思いのほか理論編が長くなってしまいすでにインパルス応答で4記事目ですが、ようやく実装の話です。
グレンジャー因果性と同様に、statsmodelsで学習したVARモデルはインパルス応答関数のメソッドも持っています。

メソッド名は irf です。
statsmodels.tsa.vector_ar.var_model.VARResults.irf

早速やってみましょう。
モデルは、直交化インパルス応答関数の計算例の記事と同じです。
サンプルデータの生成と、パラメーターの推定まで完了しているとします。
今回は次の状態になりました。


# 学習したモデルのパラメーター
result = model.fit(maxlags=1)
print(result.summary())
"""
  Summary of Regression Results   
==================================
Model:                         VAR
Method:                        OLS
Date:           Sun, 23, Feb, 2020
Time:                     21:45:36
--------------------------------------------------------------------
No. of Equations:         2.00000    BIC:                   0.924001
Nobs:                     99.0000    HQIC:                  0.830357
Log likelihood:          -312.903    FPE:                    2.15278
AIC:                     0.766721    Det(Omega_mle):         2.02800
--------------------------------------------------------------------
Results for equation y1
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         0.177695         0.918241            0.194           0.847
L1.y1         0.553623         0.112916            4.903           0.000
L1.y2         0.169688         0.169945            0.998           0.318
========================================================================

Results for equation y2
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         1.306800         0.408504            3.199           0.001
L1.y1         0.085109         0.050234            1.694           0.090
L1.y2         0.760450         0.075605           10.058           0.000
========================================================================

Correlation matrix of residuals
            y1        y2
y1    1.000000  0.629788
y2    0.629788  1.000000
"""

サマリーの最後に、撹乱項の相関行列がありますが、分散共分散行列がみたいのでそちらも表示しておきましょう。
sigma_uに入ってます。


# 分散共分散行列
print(result.sigma_u)
"""
[[4.12100608 1.1546159 ]
 [1.1546159  0.81561201]]
"""

さて、まずは非直行化インパルス応答関数を計算して可視化してみましょう。
ちなみに、直行化/非直行化を指定するのは可視化(plot)の引数であり、計算時点で両方とも算出されています。


import matplotlib.pyplot as plt

# 10期先までの非直交化インパルス応答関数のプロット
period = 10
irf = result.irf(period)
irf.plot()
plt.show()

結果がこちらです。

想定通りの結果になっていますね。

続いて、(直行化)インパルス応答関数です。
用語としては、単にインパルス応答関数と言ったら、直行化の方をさすのですが、ライブラリでは明示的に、
orth=Trueを指定しなければなりません。(デフォルトはFalse)

計算は完了しているので、可視化のみです。


irf.plot(orth=True)
plt.show()

結果がこちら。

0期先の値が特徴的で、 y2からy1へは影響していないのに、y1からy2へは影響が発生していますね。
また、もう一点注意しないといけない点があります。
非直行化の方は、y1からy1 および y2からy2への0期先の影響の値が1であることから分かる通り、
1単位のショックを与えた時のインパルス応答関数を計算されていました。
しかし、直行化の方は、そうではなく、1標準偏差のショックを与えた時のインパルス応答関数の値になっています。
(この辺の挙動はドキュメントのどこに書いてあるのだろう?)

結果をプロットするのではなく、値を利用したい場合もあると思います。
その場合、次の2変数にプロットされた値がそれぞれ入っていますのでこれが使えます。


irf.irfs
irf.orth_irfs

直交化インパルス応答関数の計算例

今回は直交化インパルス応答関数を実際に計算してみます。

例として扱うのは、非直交化インパルス応答関数の記事と同じ次の例です。
$$
\left\{\begin{matrix}\
y_{1t}&=& -1+ 0.6 y_{1,t-1}+ 0.3 y_{2,t-1}+\varepsilon_{1t}\\\
y_{2t}&=& 1 + 0.1 y_{1,t-1}+ 0.8 y_{2,t-1}+\varepsilon_{2t}\
\end{matrix}\
\right.\
,\left(\begin{matrix}\varepsilon_{1t}\\\varepsilon_{2t}\end{matrix}\right)\sim W.N.(\boldsymbol{\Sigma})
$$
$$
\boldsymbol{\Sigma} = \left(\begin{matrix}4&1.2\\1.2&1\end{matrix}\right)
$$

さて、まずは$\boldsymbol{\Sigma}$を分解する必要があります。

$$
\begin{align}
\mathbf{A}=
\left(\begin{matrix}1&0\\0.3&1\end{matrix}\right),\\
\mathbf{D}= Var(\mathbf{u}_t)
\left(\begin{matrix}4&0\\0&0.64\end{matrix}\right)
\end{align}
$$
とすると、
$$
\boldsymbol{\Sigma} = \mathbf{ADA^{\top}}
$$
が成立します。
このため、$\boldsymbol{\varepsilon_t}$は次のように、
直交化撹乱項に分解できます。
$$
\left\{\begin{align}\
\varepsilon_{1t}&=u_{1t}\\\
\varepsilon_{2t}&=0.3u_{1t}+u_{2t}\
\end{align}\right.
$$
この撹乱項を用いることで、元のモデルがこのように変わります。

$$
\left\{\begin{matrix}\
y_{1t}&=& -1&+ &0.6 y_{1,t-1}&+ &0.3 y_{2,t-1}&+&u_{1t}\\\
y_{2t}&=& 1 &+ &0.1 y_{1,t-1}&+ &0.8 y_{2,t-1}&+&0.3u_{1t}+u_{2t}\
\end{matrix}\
\right.\
$$
$$
\mathbf{u}_t \sim W.N.(\mathbf{D})
$$
$$
\mathbf{D} = \left(\begin{matrix}4&0\\0&0.64\end{matrix}\right)
$$

この形に変形できれば、非直交化インパルス応答関数と同じ様に計算が出来ます。
例えば、$y_1$に1単位のショックを与えたときの同時点におけるインパルス応答は次になります。
$$
\begin{align}
IRF_{11}(0) = \frac{\partial y_{1t}}{\partial u_{1t}} = 1\\
IRF_{21}(0) = \frac{\partial y_{2t}}{\partial u_{1t}} = 0.3
\end{align}
$$
$IRF_{21}(0)\neq0$となるのが大きな特徴です。

ただし、逆に$y_2$に1単位のショックを与えても、$y_1$には影響がありません。
要するに次のようになります。
$$
\begin{align}
IRF_{12}(0) = \frac{\partial y_{1t}}{\partial u_{2t}} = 0\\
IRF_{22}(0) = \frac{\partial y_{2t}}{\partial u_{2t}} = 1
\end{align}
$$
これが三角分解の仮定により発生している現象であり、変数の並べ方が結果に影響してしまっている点です。
このためインパルス応答関数を使う時は、変数の並べ方に気を使う必要があります。

1期後先以降の値は、非直交化インパルス応答関数のときと同じ漸化式で逐次的に求めることが可能です。

直交化インパルス応答関数

前回の記事: 非直交化インパルス応答関数
の続きです。今回も引用元は沖本本です。

最後に書いたとおり、非直交化インパルス応答関数には、撹乱項の間の相関を無視しているという問題がありました。
これを解決するには、撹乱項を互いに無相関な撹乱項に分解して、それぞれにショックを与えたときの各変数の変化を調べることです。
しかし、一般に撹乱項を互いに無相関な撹乱項に分解することは不可能で、何かしらの仮定が必要になります。
そこで、撹乱項の分散共分散行列$\boldsymbol{\Sigma}$の三角分解を用いて、撹乱項を互いに無相関な撹乱項に分解できると仮定します。

定義 (直交化インパルス応答関数):
撹乱項の分散共分散行列の三角分解を用いて、撹乱項を互いに無相関な撹乱項に分解したとき、
互いに無相関な撹乱項は直交化撹乱項といわれる。また、$y_j$の直交化撹乱項に1単位または1標準偏差のショックを与えた時の、
$y_i$の直交化インパルス応答を時間の関数としてみたものは、$y_j$に対する$y_i$の直交化インパルス応答関数といわれる。

一般的に、インパルス応答(関数)というというと直交化インパルス応答(関数)を指すそうです。
(後日の記事でstatsmodelsによる実装紹介しますが、このモジュールでは非直交化のほうをデフォルトにしてるように見えます。
もしそうだとしたら気をつけて使わないといけないので、記事にする前にもう一回検証します。)

さて、もう少し具体的に定義します。
まず撹乱項 $\boldsymbol{\varepsilon}_t$ の分散共分散行列$\boldsymbol{\Sigma}$は正定値行列であるので、
$\mathbf{A}$を対角成分が1に等しい下三角行列、$\mathbf{D}$を対角行列として、
$$
\boldsymbol{\Sigma} = \mathbf{ADA^{\top}}
$$
と三角分解することが出来ます。
(と、本には書いてありますが、分散共分散行列が正定値なだけでなく、対称行列であることも考慮する必要あると思います。)

このとき、直交化撹乱項 $\mathbf{u}_t$は、
$$
\mathbf{u}_t = \mathbf{A}^{-1}\boldsymbol{\varepsilon}_t
$$
で定義され、各変数固有の変動を表すとみなされます。
なお、
$$
\begin{align}
Var(\mathbf{u}_t) &= Var(\mathbf{A}^{-1}\boldsymbol{\varepsilon}_t) = \mathbf{A}^{-1} Var(\boldsymbol{\varepsilon}_t)(\mathbf{A}^{-1})^{\top} \\
& =\mathbf{A}^{-1}\boldsymbol{\Sigma}(\mathbf{A}^{\top})^{-1}=\mathbf{D}
\end{align}
$$
であることから、$\mathbf{u}_t$が互いに無相関であることが分かります。
また、それぞれの分散も確認できます。

$y_j$のショックに対する$y_i$の(直交化)インパルス応答関数は、$u_{jt}$に1単位のショックを与えたときの、
$y_i$の反応を時間の関数と見たものであり、
$$
IRF_{ij}(k) = \frac{\partial y_{i,t+k}}{\partial u_{jt}}, \ \ \ k=0, 1, 2, \dots.
$$
と表されます。
1単位ではなくて、1標準偏差のショックを与えたときの変数の反応を見たい時は、
$IRF_{ij}(k)$に$\sqrt{Var(u_{jt})}$を掛けることで求まります。(非直交化のときと同じですね)

1標準偏差のショックを与えたときのインパルス応答関数を求める方法として、
三角分解の変わりにコレスキー分解を用いて撹乱項を分解して、
その分解した撹乱項に1単位ショックを与える方法も紹介されていますが、
撹乱項の絶対値が標準偏差倍違うだけの話なので省略します。

非直交化インパルス応答関数

グレンジャー因果性で複数の過程の間の関係を調べる方法を紹介してきましたが、グレンジャー因果性には
関係の強さを測れないという問題がありました。
そこで、それを補うツールとしてインパルス応答関数(IRF、 impulse response function)というものがあります。

例によって、沖本先生の本を参照しながら紹介させていただきます。

インパルス応答関数はある変数にショックを与える(要するに少し値を変動させる)と、
それが他の変数に与える影響を分析することができるものです。
VARモデルにおいては、ショックの識別の仕方によって、複数のインパルス応答関数が存在しますが、
まずは非直交化インパルス応答関数について説明します。

一般的な$n$変量のVAR(p)モデルを考えます。
$$
\mathbf{y}_t = \mathbf{c}+\boldsymbol{\Phi}_1\mathbf{y}_{t-1}+\cdots+\boldsymbol{\Phi}_p\mathbf{y}_{t-p}+\boldsymbol{\varepsilon}_t,
\ \ \ \boldsymbol{\varepsilon}_t\sim W.N.(\boldsymbol{\Sigma})
$$
本では、$\boldsymbol{\Sigma}$は対角行列ではないと仮定されていますが、対角行列の場合も同じだと思います。
(対角行列だと、その直後に登場する直行化インパルス応答関数と全く同じになってしまうで説明の都合上仮定されているのかと。)

この時、非直交化インパルス応答関数は次のように定義されます。

定義 (非直交化インパルス応答関数):
$y_{jt}$の撹乱項$\varepsilon$だけに、1単位(または、1標準偏差)のショックを与えた時の、
$y_{i,t+k}$の値の変化は、$y_j$のショックに対する$y_i$の$k$期後の非直行化インパルス応答と呼ばれる。
また、それを$k$の関数としてみたものは、$y_j$のショックに対する$y_i$の非直行化インパルス応答関数と言われる。

変化量なので、偏微分を用いて次のように表記できます。(1単位のショックを与えた場合。)
$$
IRF_{ij}(k) = \frac{\partial y_{i,t+k}}{\partial \varepsilon_{jt}}, \ \ \ k=0, 1, 2, \dots.
$$
1標準偏差のショックを与えた場合は、$IRF_{ij}(k)$に$\sqrt{Var(\varepsilon_{jt})}$を掛けることで求まります。

さて、具体的な計算方法ですが、$k=0$から逐次的に計算していくことで簡単に求まります。

次のモデルを例に、計算例が乗っているのでやってみましょう。
$$
\left\{\begin{matrix}\
y_{1t}&=& -1+ 0.6 y_{1,t-1}+ 0.3 y_{2,t-1}+\varepsilon_{1t}\\\
y_{2t}&=& 1 + 0.1 y_{1,t-1}+ 0.8 y_{2,t-1}+\varepsilon_{2t}\
\end{matrix}\
\right.\
,\left(\begin{matrix}\varepsilon_{1t}\\\varepsilon_{2t}\end{matrix}\right)\sim W.N.(\Sigma)
$$
$$
\Sigma = \left(\begin{matrix}4&1.2\\1.2&1\end{matrix}\right)
$$

$y_1$に1単位のショックを与えた時の非直交化インパルス応答関数を具体的に計算します。
k=0の場合は簡単で、そのまま微分するだけです。
$$
\begin{align}
IRF_{11}(0) = \frac{\partial y_{1t}}{\partial \varepsilon_{1t}} = 1\\
IRF_{21}(0) = \frac{\partial y_{2t}}{\partial \varepsilon_{1t}} = 0
\end{align}
$$
この後は、多変数関数の合成関数の連鎖律を使って計算していきます。$k=1$の場合、$y_1$,$y_2$それぞれについて計算すると次のようになります。
$$
\begin{align}
IRF_{11}(1) &= \frac{\partial y_{1,t+1}}{\partial \varepsilon_{1t}}
= 0.6\times \frac{\partial y_{1t}}{\partial \varepsilon_{1t}} + 0.3 \times \frac{\partial y_{2t}}{\partial \varepsilon_{1t}}\\
&= 0.6 IRF_{11}(0) + 0.3 IRF_{21}(0) =0.6
\end{align}
$$
$$
\begin{align}
IRF_{21}(1) &= \frac{\partial y_{2,t+1}}{\partial \varepsilon_{1t}}
= 0.1\times \frac{\partial y_{1t}}{\partial \varepsilon_{1t}} + 0.8 \times \frac{\partial y_{2t}}{\partial \varepsilon_{1t}}\\
&= 0.1 IRF_{11}(0) + 0.8 IRF_{21}(0) =0.1
\end{align}
$$

この後の $k=1,2,3,…$についても次の漸化式で逐次的に計算できます。
$$
\left\{\
\begin{align}\
IRF_{11}(k) &= 0.6IRF_{11}(k-1) + 0.3IRF_{21}(k-1)\\\
IRF_{21}(k) &= 0.1IRF_{11}(k-1) + 0.8IRF_{21}(k-1)\
\end{align}\
\right.\
$$
この式の中に元のモデルの式の定数項ができませんが、それはインパルス応答関数が変化量に着目しているため、
無関係(微分で消えてしまった)からです。

さて、これで非直交化インパルス応答関数が計算できましたが、
これには一つ問題点があります。
というのも、最初のモデルの式の$\boldsymbol{\Sigma}$を見るとわかるとおり、$\varepsilon_{1t}$と$\varepsilon_{2t}$には
相関があることが仮定されています。
分散共分散行列から相関係数を計算すると
$$
Corr(\varepsilon_{1t}, \varepsilon_{2t}) = \frac{1.2}{\sqrt(4\times 1)} = 0.6
$$
となり、そこそこ強い相関です。
そのため、$\varepsilon_{1t}$にだけショックを与えて、$\varepsilon_{2t}$はそのままという、
非直交化インパルス応答関数には少し不自然な点が残ります。
すでに長くなってきたので、この点を解消する方法につては次の記事で紹介させていただきます。

statsmodels でグレンジャー因果性の確認

前回の記事でグレンジャー因果性という用語を紹介したので、今回はPythonで実際に計算してみましょう。
参考: グレンジャー因果性

データは、VARモデルのパラメーター推定を行った時のように、乱数で生成します。
参考: Python(statsmodels)でVARモデルのパラメーター推定

コードの重複する部分は省略して、変数resultに、推定したモデルが格納されているとします。
(今回の記事用に実行し直したので、パラメーターは前回と違います。)

この推定結果のオブジェクトが、 test_causality というメソッドを持っており、
それを使って、グレンジャー因果性の検証ができます。
ドキュメント: statsmodels.tsa.vector_ar.var_model.VARResults.test_causality


statsmodelsにはこの他にも、二つの過程のグレンジャー因果性を確認するツールとして、
grangercausalitytests というのが用意されているようなのですが、今回はこちらは使いません。

それでは早速やってみましょう。


# 推定結果のVARモデルの確認。(以前の記事のコードで生成されている前提)
print(result.summary())
"""
  Summary of Regression Results   
==================================
Model:                         VAR
Method:                        OLS
Date:           Sun, 09, Feb, 2020
Time:                     16:37:19
--------------------------------------------------------------------
No. of Equations:         2.00000    BIC:                    1.41978
Nobs:                     99.0000    HQIC:                   1.32614
Log likelihood:          -337.444    FPE:                    3.53438
AIC:                      1.26250    Det(Omega_mle):         3.32953
--------------------------------------------------------------------
Results for equation y1
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         1.263791         0.705621            1.791           0.073
L1.y1         0.604660         0.124431            4.859           0.000
L1.y2         0.267836         0.145104            1.846           0.065
========================================================================

Results for equation y2
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const        -3.289405         0.292191          -11.258           0.000
L1.y1         0.601292         0.051526           11.670           0.000
L1.y2         0.263057         0.060086            4.378           0.000
========================================================================

Correlation matrix of residuals
            y1        y2
y1    1.000000  0.596470
y2    0.596470  1.000000
"""

# y2 から y1 へのグレンジャー因果性の確認
test_result = result.test_causality(caused="y1", causing="y2")
print(test_result.summary())
"""
Granger causality F-test. H_0: y2 does not Granger-cause y1. Conclusion: fail to reject H_0 at 5% significance level.
==============================================
Test statistic Critical value p-value    df
----------------------------------------------
         3.407          3.890   0.066 (1, 192)
----------------------------------------------
"""

# y1 から y2 へのグレンジャー因果性の確認
test_result = result.test_causality(caused="y2", causing="y1")
print(test_result.summary())
"""
Granger causality F-test. H_0: y1 does not Granger-cause y2. Conclusion: reject H_0 at 5% significance level.
==============================================
Test statistic Critical value p-value    df
----------------------------------------------
         136.2          3.890   0.000 (1, 192)
----------------------------------------------
"""

当然データによって結果が変わる(乱数で生成しているので、モデルが同じでも実行するたびに変わります。)のですが、
y1からy2へのグレンジャー因果性がないという過程は棄却されているので、y1はy2に影響しているようです。
y2からy1の向きは、そうはならなかったですね。

さて、上のコード例のように、 test_causality というメソッドに、
caused と causing の引数で変数名を渡してあげるだけで、テスト結果を返してくれました。
例では変数名を指定しましたが、 0 や 1など、何番目の変数なのかを指定してもOKです。
これが2変数ではなくもっと多くの変数のVARモデルで、多変数から多変数へのグレンジャー因果性を調べる場合は、
causedやcausingをリストで指定すれば、同じようにテストできます。

さて、 test_causality は CausalityTestResults というオブジェクトを返してくれますが、
上の例のように、 summary() というメソッドを持っているので、人間の目で見る場合は、これをprintするのが良いでしょう。

p値などのプログラムで使いたい場合は、それぞれ変数に持っているのでそちらを使えます。


# 帰無仮説
print(test_result.h0)
# H_0: y1 does not Granger-cause y2

# 統計量
print(test_result.test_statistic)
# 136.18235785962648

# p値
print(test_result.pvalue)
# 3.9579973656016236e-24

# 結論
print(test_result.conclusion_str)
# Conclusion: reject H_0

グレンジャー因果性

ここ数回VARモデルの話を書いていますが、VARモデルの目的の一つは複数の時系列データ間の関連を調べることです。
普通は、各データにその背景があって、それらに関する理論からあの値がこうなったらどういう影響があってこっちの値がどうなる、と調べるのですが、
データだけから因果性の有無を判別することができると非常に便利です。
そのような考えを元に、 Grangerが提案したのがグレンジャー因果性です。

今回の記事も、沖本先生の経済ファイナンスデータの軽量時系列分析を元に各定義を紹介していきます。

定義(グレンジャー因果性)
現在と過去の$x$の値だけに基づいた将来の$x$の予測と、現在と過去の$x$と$y$の値に基づいた将来の$x$の予測を比較して、
後者のMSEの方が小さくなる場合、$y_t$から$x_t$へのグレンジャー因果性(Granger causality)が存在すると言われる。

これは、$x,y$の2つの過程に対する定義ですが、もっと多くの変数が存在する場合に拡張できます。
それが次の一般的なグレンジャー因果性です。

定義(一般的なグレンジャー因果性)
$\mathbf{x}_t$と$\mathbf{y}_t$をベクトル過程とする。また、$\mathbf{x}$と$\mathbf{y}$の現在と過去の値を含む、時点$t$において利用可能な情報の集合を
$\boldsymbol{\Omega}_t$とし、$\boldsymbol{\Omega}_t$から現在と過去の$\mathbf{y}$を取り除いたものを、$\tilde{\boldsymbol{\Omega}}_t$とする。
この時、$\tilde{\boldsymbol{\Omega}}_t$に基づいた将来の$\mathbf{x}$の予測と、$\boldsymbol{\Omega}_t$に基づいた将来の$\mathbf{x}$の予測を比較して、
後者のMSEのほうが小さくなる場合、$\mathbf{y}_t$から$\mathbf{x}_t$へのグレンジャー因果性が存在すると言われる。
ここで、MSEの大小は行列の意味での大小であることに注意されたい。

沖本本によると、グレンジャー因果性が頻繁に用いられるようになったたった一つの理由は、VARの枠組みでは、
グレンジャー因果性の有無が比較的簡単に明確な形で検証できるからだそうです。

2変量VAR(2)モデルを用いた説明が紹介されているのでみていきましょう。
モデルを具体的に書くと次のようになります。
$$
\left\{\begin{matrix}\
y_{1t}=c_1+\phi_{11}^{(1)}y_{1,t-1}+\phi_{12}^{(1)}y_{2,t-1}+\phi_{11}^{(2)}y_{1,t-2}+\phi_{12}^{(2)}y_{2,t-2}+\varepsilon_{1t},\\\
y_{2t}=c_2+\phi_{21}^{(1)}y_{1,t-1}+\phi_{22}^{(1)}y_{2,t-1}+\phi_{21}^{(2)}y_{1,t-2}+\phi_{22}^{(2)}y_{2,t-2}+\varepsilon_{2t}\
\end{matrix}\
\right.\
$$

ここで、$y_{2t}$から、$y_{1t}$へのグレンジャーレンジャー因果性が”存在しない”ということjは、
$\phi_{12}^{(1)}=\phi_{12}^{(2)}y_{2,t-2}=0$ということ同値です。
一般のVAR(p)において、$y_{2t}$から、$y_{1t}$へのグレンジャー因果性が存在しないとは、VARの$y_1$の式において、$y_2$に関係する係数が全部$0$ということになります。
その結果、VARの枠組みだと、F検定を用いてグレンジャー因果性を検定できます。
(複数の回帰係数をまとめて検定するので、t検定ではなくF検定だそうです。この辺り、自分も理解が浅いので別途勉強し直そうと思います。)

手順としては次のようになります。
$H_0: \phi_{12}^{(1)}=\phi_{12}^{(2)}y_{2,t-2}=0$ を検定すれば良いので、まずは、
$$
y_{1t}=c_1+\phi_{11}^{(1)}y_{1,t-1}+\phi_{12}^{(1)}y_{2,t-1}+\phi_{11}^{(2)}y_{1,t-2}+\phi_{12}^{(2)}y_{2,t-2}+\varepsilon_{1t}
$$
をOLSで推定して、残渣平方和を$SSR_1$とします。
次に制約を課したモデル、
$$
y_{1t}=c_1+\phi_{11}^{(1)}y_{1,t-1}+\phi_{11}^{(2)}y_{1,t-2}+\varepsilon_{1t}
$$
をOLSで推定して、その残差平方和を$SSR_0$とします。
この時、F統計量が次の式で定義されます。
$$
F\equiv \frac{(SSR_0-SSR_1)/2}{SSR_1/(T-5)}
$$
そうすると、$2F$は漸近的に$\chi^2(2)$に従うことが知られています。
なので、$2F$の値を$\chi^2(2)$の$95%$店と比較して、$2F$の方が大きければ、$y_2$から$y_1$へのグレンジャー因果性が存在しないという
帰無仮説を棄却して、$y_2$は$y_1$の将来を予測するのに有用であると結論づけられます。

これは1つの過程から別の1つの過程への影響を調べるものですが、もっと一般には次の手順になります。
(ここも沖本先生の本の丸写し)

n変量VAR(p)モデルにおけるグレンジャー因果性検定の手順
(1) VARモデルにおける$y_{kt}$のモデルをOLSで推定し、その残差平方和を$SSR_1$とする。
(2) VARモデルにおける$y_{kt}$のモデルに制約を課したモデルをOLSで推定し、その残差平方和を$SSR_0$とする。
(3) F統計量を
$$
F\equiv \frac{(SSR_0-SSR_1)/r}{SSR_1/(T-np-1)}
$$
で計算する。ここで、$r$はグレンジャー因果性検定に必要な制約の数である。
(4) $rF$を$\chi^2(r)$の$95%$点と比較し、$rF$のほうが大きければ、ある変数(群)から$y_{kt}$へのグレンジャー因果性は存在し、小さければグレンジャー因果性は存在しないと結論する。

あら?帰無仮説が棄却できなかった時に、存在しないって結論していいのかな?ここの記述はちょっと怪しいですね。

最後に、グレンジャー因果性の長所と短所がまとめられています。(沖本本で株式市場の例が紹介された後。)
長所は定義が明確であることt、データから容易に検定できることです。

短所は、まずはグレンジャー因果性が通常の因果性と異なっていることです。
グレンジャー因果性は因果性が存在する必要条件ではありますが、十分条件でありません。
なので、グレンジャー因果性があったからといって因果性があるとは結論づけられません。
また、グレンジャー因果性の方向と、通常の因果性の方向が同じになると限らず、
極端な例では逆向きだけ存在することもあるそうです。
そのため、何か別の理論等により因果関係の方向がはっきりしている時以外は、
定義通りに、予測に有用かどうかだけの観点で解釈することが良さそうです。

また、グレンジャー因果性は訂正的概念で、関係の強さは測れません。
(F統計量が大きいほど関係が強いみたいな考えは誤りということですね。)

Python(statsmodels)でVARモデルのパラメーター推定

今回は2変量のVARモデルにしたがうデータから、Pythonを使ってモデルの係数を推定する方法を紹介します。
答え合わせをしたいので、使うデータは特定のVARモデルから生成します。

モデルはいつもの沖本先生の本のP102ページの演習問題から拝借しました。

$$
\left\{\begin{matrix}\
y_{1t} &=& 2+ 0.5 y_{1,t-1}+ 0.4 y_{2,t-1}+\varepsilon_{1t},\\\
y_{2t} &=& -3 + 0.6 y_{1,t-1}+ 0.3 y_{2,t-1}+\varepsilon_{2t}\
\end{matrix}\
\right.\
\left(\begin{matrix}\varepsilon_{1t}\\\varepsilon_{2t}\end{matrix}\right)\sim W.N.(\mathbf{\Sigma})
$$
$$
\Sigma = \left(\begin{matrix}4&1.2\\1.2&1\end{matrix}\right)
$$

まずはデータを生成します。


import numpy as np
import matplotlib.pyplot as plt

mean = np.array([0, 0])
cov = np.array([[4, 1.2], [1.2, 1]])

# 少し長めの系列でデータを生成する。
data = np.zeros([110, 2])
epsilons = np.random.multivariate_normal(mean, cov, size=110)
print(np.cov(epsilons, rowvar=False))
"""
[[3.96508657 1.1917421 ]
 [1.1917421  1.01758275]]
"""

for i in range(1, 110):
    data[i, 0] = 2 + 0.5*data[i-1, 0] + 0.4*data[i-1, 1] + epsilons[i, 0]
    data[i, 1] = -3 + 0.6*data[i-1, 0] + 0.3*data[i-1, 1] + epsilons[i, 1]

# 初めの10項を切り捨てる
data = data[10:]
# data.shape == (100, 2)

fig = plt.figure(figsize=(12, 5), facecolor="w")
ax = fig.add_subplot(111)
ax.plot(data[:, 0], label="$y_1$")
ax.plot(data[:, 1], label="$y_2$")
ax.set_title("サンプルデータ")
plt.legend()
plt.show()

生成されたデータをプロットしたのがこちら。

さて、いよいよ statsmodelsでVARモデルのパラメーター推定をやりますが、
方法は非常に簡単で、時系列モデルのapiからVARモデルを呼び出してデータを渡し、学習するだけです。

ドキュメント:
Vector Autoregressions tsa.vector_ar


from statsmodels.tsa.api import VAR
model = VAR(data)
result = model.fit(maxlags=1)
print(result.summary())
"""
  Summary of Regression Results
==================================
Model:                         VAR
Method:                        OLS
Date:           Mon, 03, Feb, 2020
Time:                     01:10:12
--------------------------------------------------------------------
No. of Equations:         2.00000    BIC:                    1.23044
Nobs:                     99.0000    HQIC:                   1.13680
Log likelihood:          -328.071    FPE:                    2.92472
AIC:                      1.07316    Det(Omega_mle):         2.75521
--------------------------------------------------------------------
Results for equation y1
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         1.879511         0.616917            3.047           0.002
L1.y1         0.464346         0.124135            3.741           0.000
L1.y2         0.462511         0.132200            3.499           0.000
========================================================================

Results for equation y2
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const        -2.683660         0.305923           -8.772           0.000
L1.y1         0.561468         0.061557            9.121           0.000
L1.y2         0.385671         0.065556            5.883           0.000
========================================================================

Correlation matrix of residuals
            y1        y2
y1    1.000000  0.583971
y2    0.583971  1.000000
"""

今回はモデルの形が事前にわかっているので、 maxlags=1 としましたが、 ここを少し大きめの値にしておけば自動的にlagの値を選んでくれます。
選び方は 引数”ic”で指定できます。 (未指定ならAIC)

ic : {‘aic’, ‘fpe’, ‘hqic’, ‘bic’, None}
Information criterion to use for VAR order selection.
aic : Akaike
fpe : Final prediction error
hqic : Hannan-Quinn
bic : Bayesian a.k.a. Schwarz

さて、結果の係数ですが、summaryを見ると大体次の値になっています。(四捨五入しました)
$$
\left\{\begin{matrix}\
y_{1t} &=& 1.89+ 0.46 y_{1,t-1}+ 0.46 y_{2,t-1}\\\
y_{2t} &=& -2.68 + 0.56 y_{1,t-1}+ 0.39 y_{2,t-1}\
\end{matrix}\
\right.\
$$
とても正確というわけではありませんが、概ね正しく推定できていそうです。

ベクトルホワイトノイズ

前回の記事で弱定常過程の定義をベクトルに拡張しましたが、今回はホワイトノイズをベクトルに拡張した
ベクトルホワイトノイズを定義します。
参考: 定常過程の例としてのiid過程とホワイトノイズ

ベクトル過程 $\boldsymbol{\varepsilon}_t$ が全ての時点 $t$ において次の二つの式を満たすとします。
$$
\begin{align}
E(\boldsymbol{\varepsilon}_t) &= \mathbf{0}&\\
E(\boldsymbol{\varepsilon}_t \boldsymbol{\varepsilon}_{t-k}^{\top}) &=
\left\{
\begin{matrix}
\boldsymbol{\Sigma},& k=0\\\
\mathbf{0},& k\neq 0 \
\end{matrix}
\right.
\end{align}
$$

この時、 $\boldsymbol{\varepsilon}_t$ をベクトルホワイトノイズと呼びます。
このベクトルホワイトノイズは弱定常であり、自己相関を持たないことはすぐわかります。
$k\neq 0$の時、$k$次の自己共分散行列が$\mathbf{0}$ですから。

沖本先生の本で注意されている通り、同時点での各変数は相関を持つことができます。

ただ、一点、次の記載があるのには少し引っかかります。(初版第13刷 P.76)

$\boldsymbol{\Sigma}$ は$n\times n$正定値行列であり、対角行列である必要はないことには注意が必要である.

$n\times n$正定値行列であることは確かなのでしが、これ、必ず対角行列になるような気がします。
というのも、$\boldsymbol{\Gamma}_k = \boldsymbol{\Gamma}^{\top}_{-k}$に、$k=0$ を代入したら、
$\boldsymbol{\Gamma}_0 = \boldsymbol{\Gamma}^{\top}_{0}$ となり、これは対角行列ですよね。

僕が何か勘違いしているのかもしれませんがおそらく誤植か何かのように思います。
(出版社サイトの今日時点の正誤表には入ってないようです)