最近、カルマンフィルタにはまっている。適当なモデルを作ってデータを放り込むと、ごにょごにょして、尤もらしい答えを出してくれるところが最高である。はじめは、自分ででっちあげた疑似問題を解いていたのだが、そろそろリアルデータで試してみたいなあ、と思うようになり、表題の通り、GDPのデータを分析してみることにした。
トレンド推定
トレンドのモデルには色々あるが、今回は下記の状態空間モデルを用いることにした。
引用元：https://www.statsmodels.org/dev/examples/notebooks/generated/statespace_local_linear_trend.html
何やら難しそうに見えるが、トレンドのレベルと変化率をそれぞれμ、νとして、νはランダムウォーク、μはνの積算値＋ノイズのようにしてトレンドをモデリングしている（教科書には二階階差のトレンドモデルとほぼ同じ、みたいなことが書かれていた。計算するとそうなるっぽい）。速度と位置をxとvでモデリングするのと、式の形はほとんど同じ（だと思われる。時間幅が一定なのでdtに相当する項はないけど）。このモデルを使ったのには一応理由がある。μがGDPのレベルで、νが成長率みたいにして対応付けをするといいかも、思ったからだ。何しろ、私は商学部だったのだ。ちょっとくらい、経済知ってますよというふりをしたい。
こういうことをする場合、普通はstatsmodelsのようなライブラリを使うものだし、そうするのがいいのだが、私はいろいろ遊びたいので、自前のライブラリを作った。システムモデルをF、G、システムノイズをQ、観測行列H、観測ノイズをRとして、下記のような行列をオレオレカルマンフィルタに渡す（後からいろいろ書き替えたので、下のは疑似コード。たぶんあってると思う）。
F = np.array([[1,1],
[0,1]])
G = np.array([[0,0],[0,1]]).T
Q = np.diag([sigma1,sigma2])
H = np.array([[1,0]])
R = np.diag([1.0])
ちなみに、このGというのは上の数式にはないが、計量経済学などではこういう定式化がよく使われる。このGとQをかけ合わせて最終的なノイズ項とするのに使われるのだが、今回の場合、ζに相当するノイズはなしにしたかったのでG=[[0,0],[0,1]].Tのようにしている。このように、ノイズがのらない状態変数があるとき、このGを挟んでやると便利。というわけで、行列5個（又はGがいらないなら4個）が決まったので、これでモデルは完成。ただし、sigma1に相当するパラメータは与えてやらねばならぬ。今回はscipy.optimize.minimizeを使って、RMSEを最小化するように最適化することにした。
さて、このモデルに、米国GDPのデータを通すと、μ（GDPのレベル）だの、ν（成長率）の推計値だのが出力されるはずである。ちなみに、GDPは指数的に増加するため、前処理としてlogを取っている。では、いざカーブフィッティング。
じゃーん。
カルマンフィルタの凄い所は、こんなふうにデータを「分解」してしまえるところである。ちなみに、下の小窓に出ているのが、トレンドの速度に相当する成分で、（私の稚拙な仮定により）成長率に対応するとされる時系列成分である！
と、はじめはこれで満足だったのだが、よく見てみるとコロナのところでがっつりνの成分が変動を吸収してしまっているのが分かる。成長率なんて、一年や二年で大きく変わりそうもないもんが、こんなふうにバタバタ動くってのもねー。というわけで、工夫を試みることにした。
ペナルティをつける
さて、普通のフィルタに比べて、カルマンフィルタが優れている点は、あれこれ弄ることで「もっと、こういう条件で、フィルタして欲しいなー」みたいな改造（あるいは、改悪）ができる点である。ただし、本当なら不偏推定量を出してくれるところ、弄ったせいで統計的に偏りが生じたりもするので、そのあたりはちゃんと調べなければならない。私は、主に趣味的な動機で遊んでいるだけなのと、統計的に正しいかどうかを検証するための数学力が圧倒的に不足していることから、とりあえず、色々いじってみるだけである。
そんなわけで、「成長率」をもっと滑らかにしたいという願いを込めて、最小二乗ロスにオレオレペナルティ項をつけることにする。
def loss(x):
Yt = model(*x)
xs, Ps = do_filter(Yt, ys)
return RMSE(ys, xs[:,0]) + np.std(xs[:,1])
なんか、中途半端なコードで申し訳ないのだが、これが、scipy.optimize.minimizeに渡している関数である。minimizeはこの関数が最小になるようなx（ベクトル）を探してくれる関数。ただし、戻り値が（数値的に）微分可能でなければならない。勾配が決定できないような関数だと、上手く最適化できないので注意。ともかく、この関数の処理をざっと説明すると：
- modelでモデルを作成（xを使ってF,G,Q,H,Rを生成、カルマンフィルタYtに渡す）
- do_filterでフィルタリング＆スムージングを行う
- RMSEに、GDP Growthの標準偏差を足してロスとする
という処理を行っている。この三番目の、RMSEにGDP Growthを足すというのが、ちょっとした工夫で、ない頭を絞って考えてみたものである。で、もう一度最適化してみる。
余り変わっているようには見えない。グラフが「ちっちゃい」のでわかりにくいが、実はトレンド成分（オレンジ色の線）がちょっと滑らかになってはいる。というわけで、改善はできたが、まだアメリカの成長率は不安定なままである。このままでは米国民に申し訳ないので、もうちょっと頑張ってみる。
モデルを改善する
さて、何が悪いのかというと、私がかんがえている「せいちょうりつ」とフィルタ君が推定しようとしている値に「みすまっち」があることである。そのミスマッチは、なんとなく観測モデルが吸収してくれるような気がするが、この場合、観測ノイズの期待値がゼロなので、カルマンフィルタ君はレベル変数μとGDPの実現値の誤差を最小にしようと頑張ってしまう。それが、「みすまっち」の正体なので、それを何とかしてみることにする。
というわけで、あれこれ、小一時間考えたところで、次のような改良を思いついた。
GDPの実現値は、飽くまで実現値であり、世界には真のGDPの値が存在する。その値は、わりかしゆっくり動く、成長率という変数を足し合わせたもので決まる。一方で、実現値は、多かれ少なかれ真のGDPに近づくように変動するに違いない。
このありきたりな発想を、疑似的な数式で書くと次のようになる：
GDP_true(t) = GDP_true(t-1) + GDP_growth(t-1)
GDP_growth(t) = GDP_growth(t-1) + Noise
Economical_force(t) = alpha * (GDP_realized(t-1) - GDP_true(t-1)) + Noise(t-1)
GDP_realized(t) = GDP_realized(t-1) + Economical_force(t-1) + Noise(t-1)
GDP_observed(t) = GDP_realized(t-1) + Noise(t-1)
なんか、ずいぶんごちゃごちゃしてしまったが、こんな感じ。上の二行は先ほどのトレンドモデルと同じ。一方で、GDPの実現値は必ずしもトレンドモデルのμ（GDP_true）とは一致しない。が、もしGDPの実現値がGDPの真値からずれた場合には、それを修正するような力が働く、とかそんな「いめーじ」である。さっそく、これを「じっそう」してみることにする。
F = np.array([[0, alpha, -alpha, 0],
[1, 0, 1, 0 ],
[0, 0, 1, 1 ],
[0, 0, 0, 1 ]])
G = np.array([[1,0,0],[0,1,0],[0,0,0],[0,0,1]])
Q = np.diag([sig1, sig2, sig3])
H = np.array([[0,1,0,0]])
R = np.array([[1.0]])
カルマンフィルタは、こうやって何でもかんでも「行列表記」にしなければ動いてくれないのが玉に瑕である。さて、先ほどと同じようにオレオレペナルティ付きの損失関数で、このモデルを最適化する。
めっちゃ滑らかに成長するようになった。
当初はただ単に、偶発的なショックを吸収するようなカルマンフィルタを作りたかっただけなのだが、こうやってみてみると、結構イイ感じ！じゃない？個人的には結構満足である。実際、推定されたEconomical Forceのグラフを見ていると、意外とちゃんとアメリカのたどった軌跡が反映されている気もする。
第二次世界大戦後からぐんぐん悪くなって1970年代がどん底、1980年代のにわか景気（日本もバブルだったんだよねー、このころ。よー知らんけど）、ぱっとしない90年代、2000年代のITバブル、リーマンショック、そしてコロナとこのグラフに数十年の歴史が凝縮されているようですね（自画自賛）。成長率に関しては、ピークが二カ所あって、第二次世界大戦直後と、1970年代後半。インフレがひどかった時期かな？
とまぁ、こんな感じで、こちらの都合にあうようにデータの分解が出来ちゃうところが、カルマンフィルタの凄い所です。お笑い番組に出てくる番組プロデューサーみたいな感じで、「うーん、その成長率ってのね、フィルタちゃん、もっと、なめらかーな感じでね、うん、そのショックはドンドン・ガーっていうイメージで、何とかならない？」くらいの適当な感じでも、とりあえずの数式で表現してあげたら、それなりの結果を出してくれちゃいます。
反省点
とはいえ、こういう無茶なモデルにすると、やっぱりうまく行かないところも出てくるようで、推定されたパラメータは以下の通り：
alpha = 0.7049
sig1 = 10.0
sig2 = 8.4378
sig3 = 0.004
はっきりいって大分恣意的なので、この数値にたいした意味があるとも思えないが、私が勝手に仮定した修正力に相当する状態（Economical_Forceとかって書いたやつ）の分散が、最適化ルーチンに与えた最大値になってしまっていること。要するに、ここでつじつまの合わない部分を全部吸収しようとしている。こういうこと、やってみないとまだ分からない。ちょっと考えてみたら、こうなりそうなことくらい見当はつくと思うのだけど。
ただまぁ、目的は達せられたからこれでもいいや。実際、動作自体は思ったよりちゃんとしている気がする。今度は、運動モデルかなんかを使って、ちゃんと速度が推定されるかどうかなど試してみようかなぁ。
むすび
いやね、なにも、私のモデルが経済的に優れている、みたいなことが言いたいわけではなく、カルマンフィルタで外れ値的なものをうまく回避する方法を、ちゃんと自分で考えれるようになったよ、というアピールがしたかった。アピールポイントはそっちです。
実際、時系列データは、いまどきネットを探せばいくらでも転がっているので、分析対象には事欠かない。ただ、ARMAモデルを当てはめました、みたいなことをやったところで、それほど面白くもない。
そこでですよ！カルマンフィルタです。こんなふうに、うすぼんやりでも、ボンクラでも、ちょっと工夫するだけで、「でーたもでらー」の気分に浸れる。そんな満足をもたらしてくれるのが、カルマンフィルタ。
データサイエンティストの真似事をしたいのなら、まずはこのあたりから始めてみるのが、面白いかもしれません（そして速い。パラーメータ推定に一分もかからないもんね）
追記
あとで、もうちょっといじってみると、ちゃんとパラメータが推定できるようになった。どうやら、scipy.optimize.minmizeに渡せる種々のアルゴリズムは、モデルによっては異なる結果を返すらしい（なんで？最小化するんだから同じじゃないの？）。もしかしたら、変なペナルティ項を作ってしまったので、そのせいかも。最小化のアルゴリズムなんて、どれも同じと思っていたがそうでもないらしい。
ちなみに、これが結果。
アルゴリズムはCOBYLA。三番目のグラフに季節成分らしきものが混じっている。パラメーターは下記の通り：
array([0.50633383, 3.15376758, 2.2562785 , 0.00585066])
「とりあえず、最大値にしておけばいっか」みたいな動作ではなく、ちゃんとそれっぽい数値が出てきた。最初に試したときにはTNCというアルゴリズムを使っていたのだが、そっちの方もパラメーターの範囲に余裕を持ってやると、ちゃんとした推定値をくれる。
ただ、改めて結果をよく見てみると、こっちの方は収束しきってないっぽい。やっぱり、変なモデルを食わせてもダメか…。
>>> res
fun: 0.006028753955433918
jac: array([-9.24165254e-05, -5.65501592e-05, -5.65203219e-05, -1.52698002e-03])
message: 'Max. number of function evaluations reached'
nfev: 505
nit: 16
status: 3
success: False
x: array([6.36780243e-01, 8.37106373e+00, 8.38189911e+00, 3.31502528e-03])
>>> res.x
array([6.36780243e-01, 8.37106373e+00, 8.38189911e+00, 3.31502528e-03])
以上です。