QZ法の実装¶

講義ノートの記法を踏襲し, $x^1_t$ を先決変数, $x^2_t$ を非先決変数として, システム変数の分解 $x_t = (x^1_t, x^2_t)$ ができているとする.

システム方程式

\begin{align} Ex_{t+1} = Ax_t + Bu_t \end{align}

の解から爆発的に発散するものを取り除くことで次の方程式を得る.

\begin{align} x^1_{t+1} &= \Omega_x x_t^1 + \Omega_u u_t + \Omega_y y^u_{t+1}\\ x^2_t &= \Psi_x x_t^1 + \Psi_y y^u_t \end{align}

これが Blanchard-Kahn (1980) にはじまる合理的期待モデルの分析の味噌である.

ただし, $y^u_t = y^u(u_t, u_{t+1}, \dots)$ は未来のショックに対する期待・予見を反映して決まる.

以下では Python を使って手を動かしながら, 理論の応用方法を習得してほしい.

準備¶

サードパーティライブラリの Numpy と Scipy をインポートしておく. 行列計算のための scipy サブパッケージ scipy.linalg を利用するので, こちらもインポートしておく.

import _[LIBRARY]_ as _[ALIAS]_

という構文は頻繁に使うので, 手になじませておくべきだろう.

筆者と読者の環境の違いを把握できるようにバージョンを確認しておこう.

QZ 分解¶

数値計算のための基本となる道具は, ペンシル $(E,A)$ のQZ分解である. 同じサイズをもつ任意の正方行列 $(E, A)$ に対してユニタリ行列 $Q$, $Z$ が存在して,

\begin{align} Q^* E Z = S, \qquad Q^* A Z = T \end{align}

が上三角行列になる. 特に, $(E,A)$ がレギュラー ($\exists s\in\mathbb{C}; \det (sE-A) \neq 0$) であれば, 対角成分 $(s_{ii}, t_{ii})$ が一般固有値に関する情報をもっている. この条件がなければ均衡経路の一意性はまったく期待できない. やや技術的ではあるが, レギュラー性の仮定は以下で $T_{uu}$ が可逆であることをいうのに使っている. 各自確認しほしい.

例1) QZ分解の例¶

宿題3 [16EDH3], 問題2の行列ペア

\begin{align} E = \begin{bmatrix} 0 & -1 & 3\\ 0 & 0 & -1\\ 0 & 0 & -1 \end{bmatrix}, \qquad A = \begin{bmatrix} 2 & 2 & -2\\ -1 & 0 & 0\\ -1 & 0 & -1 \end{bmatrix} \end{align}

に対するQZ分解を計算してみよう.

まずは Numpy 配列の生成から.

Python で QZ分解を実行するには, scipy.linalg の qz あるいは ordqz を使えばよい. 一般化固有値の並べ替えが必要なので, ordqz を使うのがよいだろう. まずはドキュメントを見てみよう. ドキュメントやエラーメッセージの解読はプログラミングの大切な一部なので, 決して避けて通ろうとしないように.

ordqz には必要な情報のすべては書かれていなかった (例えば A と AA の関係がこれだけでは分からない). See also の欄にある qz のドキュメントも見ておこう.

これらの内容を合わせると, ordqz は

AA, BB, alpha, beta, Q, Z = ordqz(A, B)

の形式で呼び出し, AA が A を上三角化したもの, BB が B を上三角化したものであることが分かる. 分解は A=Q*AA*Z', B=Q*BB*Z' が成り立つように行列Q, Z が作られている. これは我々の記法とほとんど同じであるが, ordqz は実QZ分解を行なっているという点で異なっている. その事実は, Parameters の欄の

output : str {'real','complex'}, optional
    Construct the real or complex QZ decomposition for real matrices.
    Default is 'real'.

というオプション引数の説明から分かる. デフォルトは 'real' であり，あえて複素QZ分解（普通のQZ分解）を出力するには, ordqz(A, B, output='complex') の形式で呼び出す必要がある. しかし, 実QZ分解でもまったく問題がないのでこの違いは気にしないでおこう. Klein (2000) は実QZ分解を使った計算の方が高速であると述べているので, わざわざ複素QZ分解を使う理由はない.

さらに，ordqz のドキュメントには

On exit, ``(ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N``, will be the
generalized eigenvalues.

と書かれているので, このQZ分解は

\begin{align} Bx_{t+1} = Ax_t \end{align}

という線形方程式に対応しているものと分かる. したがって, われわれの記号の下では, 次のように呼び出せばよい.

divide by zero (ゼロ割) の警告は非正則な E を入れるとでてしまうが, これは安全に無視できる. 煩わしい場合は消すこともできるが, さしあたって問題にならないのでそのままにしておこう.

$T = Q^* A Z$, $S = Q^* E Z$ が成り立つことを確認しておこう. 数値計算では「完全に一致すること」を期待してはいけない. 十分近くにあるかどうかを確かめる. これは numpy.allclose() で確認できる.

ちなみに, Python では Q の共役転置行列は Q.conj().T などと書く.

これで間違いなく望む上三角化できていることが分かった.

例2) 固有値のソーティング¶

Klein の方法では左上ブロックに一般化固有値の絶対値が$1$以下のもの, 右下ブロックには絶対値が$1$より大きいものが来るようにしなければならないことを思い出そう. 一般化固有値は tii / sii で取り出せるのでこれを見ればよい.

したがって -2 という有限固有値と, 代数的重複度$2$の無限大固有値が存在することが分かる. すべてが不安定固有値なので, 安定-不安定の分解は自明なものではあるが, うまく行っている.

しかし, これはいつも成り立つ訳ではない. 次の例を見れば分かるように, ordqz のデフォルトでは固有値の並びは制御されていない. 以下では次の行列ペアを考えよう.

並び替えは上手くいっていないようだが, どのようにすれば Klein (2000) に準拠した並び替えができるだろうか？

再びドキュメントに戻って Parameters の項目を見てみると, ordqz は sort パラメータを取れることが分かる.

    Parameters
    ----------
    ...

    sort : {callable, 'lhp', 'rhp', 'iuc', 'ouc'}, optional
        Specifies whether the upper eigenvalues should be sorted.  A callable
        may be passed that, given a eigenvalue, returns a boolean denoting
        whether the eigenvalue should be sorted to the top-left (True). For
        real matrix pairs, the sort function takes three real arguments
        (alphar, alphai, beta). The eigenvalue
        ``x = (alphar + alphai*1j)/beta``.  For complex matrix pairs or
        output='complex', the sort function takes two complex arguments
        (alpha, beta). The eigenvalue ``x = (alpha/beta)``.
        Alternatively, string parameters may be used:

            - 'lhp'   Left-hand plane (x.real < 0.0)
            - 'rhp'   Right-hand plane (x.real > 0.0)
            - 'iuc'   Inside the unit circle (x*x.conjugate() < 1.0)
            - 'ouc'   Outside the unit circle (x*x.conjugate() > 1.0)

callable というのは荒く言えば関数のことで, callable をパラメータに取れるということは, 固有値の並び順を決めるソート関数を自由に定義できることを意味している.

われわれの目的では, 'iuc' という定義済み順序を使えば十分そうに見えるが, 絶対値が$1$の固有値があるときにはうまくいかないのでソート関数ioucを定義しておこう.

練習問題¶

np.ordqz(A, E, sort='iuc') の形式で QZ 分解を行うとKleinが指定した固有値の並びと合わなくなる場合がある. 例を1つ探しなさい

ソート関数 iouc¶

実際に確かめてみよう.

これで望ましい結果を得ることができた.

なお, $1$ より小さい絶対値をもつ固有値の数 $n_s$ は次のようにして取得できる.

行列ブロックの表現¶

関連する行列のブロック表現を確認しておこう.

\begin{align} Z= \begin{bmatrix} Z_{1s} & Z_{1u}\\ Z_{2s} & Z_{2u} \end{bmatrix} \qquad Q^* B = C = \begin{bmatrix} C_s \\ C_u \end{bmatrix} \end{align}\begin{align} T= \begin{bmatrix} T_{ss} & T_{su}\\ 0 & T_{uu} \end{bmatrix} \qquad S= \begin{bmatrix} S_{ss} & S_{su}\\ 0 & S_{uu} \end{bmatrix} \end{align}

引き続き先ほどの例2で使った行列ペアを考えよう. $B$ および $C$ 以外はすべて明らかになっている. 例えば, 次のように設定しよう.

さて, ブロック表現はどのように扱えばよいだろうか. まず,

$$ n = 4 $$

および

$$ n_u = n - n_s = 2 $$

であることに注意しよう. さらに,

$$ n_1 = 2, \qquad n_2 = 2 $$

であるとしておこう. すべてのブロック行列をいちいち定義する必要はなく, Numpy のスライスオブジェクトを使えばいい.

例えば $Z_{1s}$, $Z_{2u}$, $T_{ss}$, $T_{uu}$ を取得する場合には次のようにすればよい.

$y^u$ の計算¶

一般の $y^u$ は無限和

\begin{align} y_{t}^{u}=-\sum_{k=0}^{\infty}\left(T_{uu}^{-1}S_{uu}\right)^{k}T_{uu}^{-1}C_{u}u_{t+k} \end{align}

によって計算されるので, 通常, 正確な評価はできない. しかし,

• ${u_t}$ が有限列であるケース
• $u_{t+1} = \Phi u_t$ が成り立つケース

では計算可能である. 後者は, $t=0$ にショックが起きて, 徐々にショックの影響が小さくなっていく ($|\Phi|<1$) ような変動のモデルとして扱いやすい. 具体的な計算方法は講義ノートに解説を譲ろう.

ここでは前者を考える. すなわち, 有限期間でショックの影響が完全になくなるケースである. $t = 0 \dots T$ のシミュレーションを行うとして, $(u_0, \dots, u_T)$ から $(y_0^u, \dots, y_T^u)$ を計算したい. ただし, $u_{T+1}$ 以降 $u$ はゼロであるとする.

\begin{align} y_{t}^{u}=-\sum_{k=0}^{T-t}\left(T_{uu}^{-1}S_{uu}\right)^{k}T_{uu}^{-1}C_{u}u_{t+k} \end{align}

に注意すれば,

\begin{align} y_{T+1}^u = 0 \end{align}

および

\begin{align} y_t^u = T_{uu}^{-1}\left[ S_{uu} y_{t+1}^u - C_u u_t \right], \quad t=0,1,\dots ,T \end{align}

によって, $y_t^u$ の系列を決定できる.

適切なサイズのベクトル列 u (今の設定では数列) を入力に取り, $y^u$ に相当するベクトル列を吐き出す関数を定義しよう. u

係数行列の計算¶

\begin{align} \Omega_{x} &= Z_{1s}S_{ss}^{-1}T_{ss}Z_{1s}^{-1} \\ % \Omega_u &= Z_{1s}S_{ss}^{-1}\left[C_{s}-(T_{su}-T_{ss}Z_{1s}^{-1}Z_{1u})T_{uu}^{-1}C_{u}\right] \\ % \Omega_y &= (Z_{1u}-Z_{1s}S_{ss}^{-1}S_{su})+Z_{1s}S_{ss}^{-1}(T_{su}-T_{ss}Z_{1s}^{-1}Z_{1u})T_{uu}^{-1}S_{uu} \end{align}

\begin{align} \Psi_x &= Z_{2s}Z_{1s}^{-1} \\ \Psi_y &= Z_{2u}-Z_{2s}Z_{1s}^{-1}Z_{1u} \end{align}

を計算すればよい.

シミュレーション¶

このモデル $(E, A, B)$ は Hansen の RBCモデルを特定のパラメータに対して計算した結果である. 同じパラメータのもとで定常状態は次のように計算できる.

ここで作成したコードを汎用的なものにするには, 関数やクラスを作って上記の手続きを抽象化しておくとよい. 各自で試してほしい.

保存¶

使い回しのために係数行列を保存しておこう.