2021-08-29

大量のバックテストで生じる多重比較の問題~False Discovery~

統計・機械学習投資

今回はファイナンスのジャーナルThe Journal of Financeなどで、

false discovery, p-hacking, harkingなどという形で、

直近の数年くらい盛り上がっている多重比較の問題を取り上げたいと思います。

多重比較とは

多重比較について詳しいことは天下のTJO氏のブログに記載されているので

そちらをご確認いただければと思います。

tjo.hatenablog.com

もともとは医療や生物分野で使われることが多かった統計分野だと思います。

最近だとPCR検査関連の偽陽性や偽陰性なども話題になりましたね。

さて、ファイナンスの論文投稿では、独自の新ファクターを提唱したうえで、

そのファクターが資産価格の決定要因となっていると主張したりするわけですが、

その際の検証において、パラメータや検証期間などをこねくり回して、

強引に有意な結果をでっちあげてジャーナルを通すことが問題視されています。

実際、今までいろいろな人が提唱してきて有意性を示したファクターは

合計すると300やら400やらともいわれており、factor zooと揶揄されています。

ファイナンスとfalse discovery rate

ファイナンスでfalse discoveryが注目されているのは上で述べた通りなのですが、

ことファクターの文脈で顕在化している背景には、

真に有効なファクターは少ないだろうというのがポイントになってきます。

というのも、本当に陽性なサンプル(TruePositive+FalseNegative)が稀有なシチュエーションでは、

false discovery rate（FDR= $\displaystyle \frac{FalsePositive}{FalsePositive + TruePositive}$ ）

が高くなってしまう、つまり誤って有意と判断することが起こりやすくなるからです。

その点について確認していきたいと思います。

FDRと真の陽性率の関係

まず1からFDRを引いた、陽性的中率（PPV）という指標を考えます。

陽性適中率 - Wikipedia

これは、陽性（有意）反応が検出された際に、

本当に陽性である条件付き確率のことを指します。

このPPVについては、真の陽性率 $\displaystyle \pi=\frac{TruePositive+FalseNegative}{AllSample}$ を用いて以下の関係性が得られます。

$\displaystyle PPV = \frac{\pi * Se}{\pi * Se+(1-\pi) * (1-Sp)}$

$Se$ は感度を、 $Sp$ は特異度を表しています。

検査の精度を $Se = 0.7, 1-Sp = 0.01$ と固定して

(↑陽性の人をきちんと陽性と検知できる確率が0.7で、
陰性の人を誤って陽性と検知する確率が0.01（有意水準1%）ということ）

PPVと $\pi$ の関係をplotすると以下のようになります。

PPVとFDRには FDR = 1 - PPV という関係があることを思い出すと、

上図から、真の陽性率 $\pi$ が小さい領域ほどFDRが大きくなることが確認できます。

多重比較に対する対応策

多重比較に対する対処法としては、2つのアプローチがあり

FWERをコントロールする方法とFDRをコントロールする方法です。

言い換えれば、FWERをコントロールする方法とは

全体に占める偽陽性の数そのものを減らす発想で、

FDRをコントロールする方は

検査で陽性になった中での偽陽性を減らす発想ということかと思います。

FWERをコントロールする手法としてはBonferroni法が有名で、

FDRをコントロールする手法としてはBenhamini-Hochberg法が有名です。

Bonferroni 法

この方法の発想はシンプルで、

「たくさん試行することでfalse discoverしてしまうなら

試行回数に応じて閾値（有意水準）も小さくしてしまえ！」という発想です。

これをもう少し数学的に記述すると、

$i=1,2,...,K$ に対して有意水準 $\frac{\alpha}{K}$ で $H_{0,i}$ を棄却する事象を $A_i$ とし、一般性を失うことなく帰無仮説族 $\mathcal{H}_{K}$ の中の最初の $L$ 個の $H_{01},...,H_{0L}$ が正しい帰無仮説とする。このとき
P(正しい帰無仮説のうち1つ以上を棄却する事象) $= P(\bigcup_{i=1}^L{A_i}) \le \sum_{i=1}^{L}P(A_{i})=L\alpha/K \le \alpha$

となるので、 $K$ 回試行したうえで有意水準 $\alpha$ の選択をするには

各試行では有意水準 $\frac{\alpha}{K}$ でやりましょう、という手法です。

ただ、この手法は確かにFWERをバウンドできてはいるのですが厳しすぎるともいわれます。

イメージとしては、偽陽性を確かに回避できてるけど偽陰性が増えてしまって

真陽性まで潰してしまう可能性を警戒しなければいけなさそうです。

Benjamini-Hochberg 法（BH法）

BH法の手順は以下の通りです。

達成したいFDRをまず決める（ここでは $\alpha$ <1とする）
$m$ 個の帰無仮説に対して検定を行い、 $p$ 値を計算する
$p$ 値を順番に並べて $p_1$ ≤ $p_2$ ≤ ... ≤ $p_m$ とする
$j = m$ として $p_j \le \frac{j}{m}\alpha$ が満たされるかを確認する
4.が満たされていれば 6. に移行し、そうでなければ $j=j-1$ として4.に戻る
帰無仮説1,2,...,j を棄却して完了

まず、以下のリンク先の内容から、帰無仮説に基づくならその $p$ 値が一様分布に従うことが分かります。 Probability integral transform - Wikipedia
仮にすべて帰無仮説に従っていた場合、 $p_j$ は一様乱数をソートしたものなので $\frac{j}{m}$ となることが想定され、 4.の基準は満たされずに、どの仮説も棄却されないことになります。
他方で、一部が帰無仮説に従っていない場合は $p_j$ の分布が歪むことで、一部の仮説が棄却されることになります。

ちょっと言葉だけだとわかりにくいですが、

下のリンク先の図を見ると直感的に納得できるかと思います。

https://www.ieice.org/~sita/forum/article/2019/201901142020.pdf

またファイナンス文脈に戻ると

さて、上記のような多重比較に対する処方箋をカスタマイズするなりして

検証を行うと、今まで提唱されてきたファクターの大半はほぼ意味がなかったよね

というようなものが、ちょっと前のファイナンスのアカデミアで

盛り上がった内容だったのかと思います。

もっとも、冒頭のTJO氏のブログにあるようにベイジアン的に、

分布に関する豊かな情報をみて判断したほうが良いのでは、

というのはその通りかと思います。

あと、余談になりますが、ファイナンスの中でも特に投資シグナル的な観点だと

別のアプローチとしてDeflated Sharpe Ratioというアプローチもあるので

また別の機会に紹介したいと思います。

今回の関連書籍

リンク

あんまりこの分野は詳しくはないのですが、一応持ってる本のご紹介です。

2021-07-22

WIndowsでGPU版lightGBMを使えるようにするメモ

Python・R PCセットアップ

GPU（RTX 3060 ti）を積んだwindowsにて

lightGBMのGPU版を使えるようにセットアップしたので

その作業メモとなります。

参考となる公式ページはこちらです。

LightGBM/python-package at master · microsoft/LightGBM · GitHub

事前に必要なもの

Visual Studio（正確にはVS build toolsの中のC/C++パッケージ？）が必要となります。

手順

基本的には、以下で述べる3つのソフトウェア(Cmake, OpenCL, Boost)が

インストールされていればpip installして完了です。

step.1 Cmakeのインストール

cmake.org

大半の人は64bit版windowsだと思うので、

上記のページから、Windows x64 Installerを入手してください。

cmake-3.21.0-windows-x86_64.msiみたいな感じのファイル名のやつです。

それをダウンロードしたら、実行してインストールが完了です。

(winget install だとpathが設定されないのでうまくいかない可能性。

なお、すでにpath通していれば winget upgrade Kitware.Cmake で大丈夫かと)

step.2 OpenCLのインストール

公式ページを見ていくと

とあります。

NVIDIA GPUの場合は cuda toolkit を入れればよいということみたいです。

cuda toolkitの導入はこちらの記事をご参照

dw-dw-dt.hatenablog.com

step.3 Boostのインストール

以下のリンク先から、自分の環境のVisual Studioと対応した

boostのexeファイルを入手します。

(ちなみに、boostは、C/C++とかの、

標準ライブラリになる一歩手前のライブラリだったような気がします。)

Visual Studio 2015 -> msvc-14.0-64.exe
Visual Studio 2017 -> msvc-14.1-64.exe
Visual Studio 2019 -> msvc-14.2-64.exe

sourceforge.net

ダウンロードが完了したら、（警告が出るけど恐れず）実行してインストールします。

step.4 pipインストール

以上でpip installが実行可能になったかと思いますので

pip install lightgbm --install-option=--gpu

を実行すればOKです。

今回の関連図書

リンク

lightgbmといえばKaggleということで定評のあるこちらを

リンク

こちらも、もともと海外のGMの方が書かれたものの邦訳ということでおすすめです

2021-06-27

カルマンフィルタのPython実装とnumbaによる高速化のお話

Python・R 統計・機械学習

時系列系の分析をしていく中で、状態空間モデルを使うこともあるかと思います。

その際、シンプルなカルマンフィルタであってもなかなか複雑だったりして

思うようにコーディングできなかったりします。

そこで、Python実装時のforループ劇遅問題含めてまとめてみました。

（pythonでカルマンフィルタといえばPykalmanもありますが、

最近メンテされてなさそうなので...）

カルマンフィルタとは

線形・ガウス型の状態空間モデルといわれるやつで、

状態空間モデルの中でも比較的シンプルなものです。

モデルとしては、状態方程式と観測方程式の2本の式から構成されます。

$x_t = Tx_{t-1} + \epsilon_t$

$y_t = Hx_t+\eta_t$

ここで、 $x_t$ は観測できない状態変数を表しており、

観測できるのは $y_t$ だけというモデルになります。

この設定において、入手できたデータ $y_0, y_1,...,y_n$ から

$\eta, \epsilon$ の分散や $T, H$ といったパラメータ推定を行ったり、

パラメータが分かっていればその状態変数を求めたり(フィルタリング)

など、いろいろな活用法があります。

どの処理をコーディングしたか

（下記関連図書を踏襲して）予測とフィルタリングを行う関数を実装しました。

引数は

前期の状態変数のフィルタリング分布( $x_{t-1|t-1}, V_{t-1|t-1}$ )

パラメータ( $T, V[ \epsilon ], H, V[ \eta ]$ )

今期の観測変数 $y_t$

で、アウトプットは

今期の観測変数の予測値 $y_{t|t-1}$

状態変数のフィルタリング分布( $x_{t|t}, V_{t|t}$ )

となっています。

これを逐次的にデータに適用することで、

所定のパラメータ下の尤度を計算できるので、

scipy.minimizeなどで最適化することで、

最尤法でパラメータ推定などもできます。

計算ロジックの概要

上記のコーディングを行うにあたっての計算ロジックは以下の通りです。

入力

$x_{t-1|t-1}, V_{t-1|t-1}, T, H, V[ \epsilon ], V[ \eta ], y_t$

step.1 今期の状態変数の予測分布

$x_{t|t-1} = Tx_{t-1|t-1}$

$V_{t|t-1} = TV_{t-1|t-1}T^{T}+V[ \epsilon ]$

step.2 今期の観測変数の予測値

$y_{t|t-1} = H x_{t|t-1}$

$y_{t|t-1}$ の分散 $=HV_{t|t-1}H^{T}+V[ \eta ]$

step.3 今期の状態変数のフィルタリング分布

カルマンゲインという数値を以下のように設定し
$K_{t} = V_{t|t-1}H^{T}(HV_{t|t-1}H^{T}+V[ \eta ])^{-1}$

状態変数に関して以下が得られます
$x_{t|t}=x_{t|t-1}+K_{t}(y_{t}-y_{t|t-1})$

$V_{t|t}=(I-K_{t}H)V_{t|t-1}$

出力

$x_{t|t}, V_{t|t},y_{t|t-1}, HV_{t|t-1}H^{T}+V[ \eta ]$

コード

可読性が劣悪ですが、下記のようなコードを作成しました。

おそらく、この関数をfor文などでループ処理することが多いと思いますので

関数のコンパイルをよしなにやってくれるnumbaを使って高速化してあります。

関数の上に @njit を付けると早くなります。

(numbaはnumpyに特化したコンパイラです。対応してない関数もあるけど...)

import numpy as np
from numba import njit


@njit    
def Num_Dot(a: np.ndarray, b: np.ndarray):
    """
    np.dotを使えば良いけど、numbaの練習がてら
    """
    c = np.zeros((a.shape[0], b.shape[1]))
    for i in range(a.shape[0]):
        for j in range(b.shape[1]):
            for k in range(a.shape[1]):
                c[i][j] += a[i][k] * b[k][j]
    return c


@njit
def kalman_filter(
    state_variable_mean_pre: np.ndarray, state_variable_covariance_matrix_pre: np.ndarray, observation_matrix: np.ndarray,
    observation_error_covariance_matrix: np.ndarray, state_transition_matrix: np.ndarray, state_error_covariance_matrix: np.ndarray,
    observations: np.ndarray
    ):
    """
    t-1のフィルタリング分布（平均と分散共分散行列）を受け取り,t時点のフィルタリング分布（平均と分散共分散行列）を返します.
    """
    state_variable_dim = state_variable_mean_pre.shape[0]
    
    # 1期先予測状態変数
    state_variable_mean = Num_Dot(state_transition_matrix, state_variable_mean_pre)
    state_variable_covariance_matrix = Num_Dot(Num_Dot(state_transition_matrix, state_variable_covariance_matrix_pre), state_transition_matrix.T) + state_error_covariance_matrix
    
    # 1期先予測観測変数
    observation_variable_mean = Num_Dot(observation_matrix, state_variable_mean)
    observation_variable_covariance_matrix = Num_Dot(Num_Dot(observation_matrix, state_variable_covariance_matrix), observation_matrix.T) + observation_error_covariance_matrix
    
    # kalman gain
    kalman_gain = Num_Dot(Num_Dot(state_variable_covariance_matrix, observation_matrix.T), np.linalg.inv(observation_variable_covariance_matrix))
    
    # フィルタリング状態変数
    filtered_state_variable_mean = state_variable_mean + Num_Dot(kalman_gain, (observations - observation_variable_mean))
    filtered_state_variable_covariance_matrix = Num_Dot((np.eye(state_variable_dim) - Num_Dot(kalman_gain, observation_matrix)), state_variable_covariance_matrix)
    
    return filtered_state_variable_mean, filtered_state_variable_covariance_matrix, observations - observation_variable_mean, observation_variable_mean, observation_variable_covariance_matrix

今回の関連図書

リンク

Rではありますが、カルマンフィルタ部分の実装などが細かく載っており参考になりました。
またパラメタの最尤推定時のコードも載っており、method=BFGSを用いるのも参考になりました。

リンク

こちらは理論面に重点を置きながらときおりコードを交えて説明してくれています。

リンク

pythonでの状態空間モデルやHMMなどを含む時系列モデリングについて載っています。
Pykalmanについても載っています。
時系列はRのものが多いのでPythonのものは珍しい印象です。

2021-06-05

WindowsでCUDAを使えるようにするメモ~GPU版pytorch~

Python・R PCセットアップ

GPU（RTX 3060 ti）を積んだwindowsにて

pytorchを使えるようにCUDAをセットアップしたので

その作業メモとなります。

(Dockerが使えれば楽なのですが、Windows版のDockerでcuda使うためには

現時点ではinsider previewに登録しないとダメっぽいので...)

事前に必要なもの

Visual Studio（正確にはVS build toolsの中のC/C++パッケージ？）が必要となります。

手順

基本的には、以下で述べる3つのソフトウェアをインストールすればOKです。

ただし、ライブラリのバージョンなどを踏まえて

あえてバージョンを落としてインストールしたりするのが面倒ではあります。

step.1 NVIDIAドライバーの入手

www.nvidia.co.jp

上記のページから、GPUに対応するインストーラーを入手してください。

今回の3060tiの場合は以下のようになりました。

このNVIDIAドライバーはあまりバージョン対応など考えず

とりあえずデフォルトの最新版を入手で大丈夫です。

※日ごろからGeforce Experienceアプリで最新ドライバにしている人はこのプロセスは不要かもしれません。

step.2 CUDA toolkitの入手

これはまず、自分が使いたいライブラリの対応状況を確認する必要があります。

今回はpytorchですので、CUDA 11.1 が必要となりそうです。

で、先ほどのNVIDIAドライバーがCUDA 11.1に対応していることを

念のために以下リンク先で確認。

（NVIDIAドライバーは後方互換性があるみたいなので、

最新のやつなら古いCUDAにも対応しているみたいです）

docs.nvidia.com

確認が取れたら、以下のアーカイブから欲しいバージョンのCUDAを取得。

(私のバージョンは11.1.1でした)

developer.nvidia.com

step.3 cuDNNの入手

以下のリンク先からcuDNNを入手します。

進めていき、step.2で取得したCUDAのバージョンと対応するcuDNNを選べばOKです。

developer.nvidia.com

step.4 インストール

これまでで入手したものを順にインストールします。

まずはNVIDIAドライバー。

次にCUDA toolkit。ちなみにここでVisual Studioが必要になります。

Visual Studioのバージョンが違う的なメッセージが出たりしますが、

それは無視して進めて問題は出なかったです。

で、最後にcuDNNのzipファイルを展開します。

展開後は、zipの中身のcudaフォルダを適当なCドライブ直下などに配置し

ユーザー環境変数のPathにC:/cuDNN/cuda/bin通して完了です。

動作確認

以上でGPUが使えるようになったと思いますのでpytorchをインストールして動作確認です。

インストールのコマンド例は、最新版なら以下

pytorch.org

ちょっと古めのコマンド例は以下

pytorch.org

となります。私の場合、cudatoolkit=11.1なので古めのコマンドを参照し

conda install pytorch==1.8.0 torchvision==0.9.0 torchaudio==0.8.0 cudatoolkit=11.1 -c pytorch -c conda-forge

で、インストールが完了したら、

import torch
print(torch.cuda.is_available())

を実行し、エラーがなければOKです。

今回の関連図書

リンク

pytorch習熟にはこのあたりの本から始めるのがいいかなと思います。

2021-04-24

東芝SBMの量子アニーリング？でポートフォリオ最適化メモ

Python・R 投資

f:id:dw_dw_dt:20210417013820j:plain

※今回の記事は専門じゃないことを記載しているので、間違いがあるかもですが、その際はコメントいただけると幸いです。

さて、今更感がありますが、東芝のシミュレーテッド分岐マシン(以下SBM)なるソルバーを使ってみました。

www.global.toshiba

このSBMというのは量子コンピュータ（イジング方式）が得意な計算を

普通のコンピュータでできるようにしてみたアルゴリズムのことだそうです。

イジング方式は組み合わせ最適化問題が得意（超早い）ということで

実用性はいったん置いておいて、とりあえずポートフォリオ最適化を

整数計画にアレンジした際のやり方をメモしました。

ノリ的にはalpacaさんがarXivにあげている以下の論文と似ているかと思います。

https://arxiv.org/pdf/2009.08412.pdf

なお、SBMは以下のAWS marketplacesにて販売されています(FPGAではなくGPU版)。

aws.amazon.com

量子アニーリング

余裕で専門外なのでざっくりとした話しかできないのですが

そもそも量子コンピュータというのは2種類に分類できるもので

ゲート方式

イジング方式

という2種類があるみたいです。

で、ゲート方式というのは汎用的なマシンのことで

イジング方式は組み合わせ最適化に特化した量子コンピュータのようです。

IBM-Qやgoogleマシンがゲート方式で、D-WAVEがイジング方式みたいです。

で、そのイジング方式は、量子アニーリングという手法で

量子焼きなまし法 - Wikipedia

組み合わせ最適化を解くのですが、それを普通のPCで解けるようにしたのが

今回触ってみる東芝様のSBMというものになります。

ちなみに量子コンピュータそのものについては、

東大が公開している以下のmaterialが役に立ちそうです。

utokyo-icepp.github.io

問題設定

あまり筋が良いとは思えませんが、

いわゆるポートフォリオ最適化に対してSBMを適用してみたいと思います。

通常のポートフォリオ最適化は

$w \in \mathbb{R}$ を最適ウェイトとして $max_{w} \mu^{T} w - \lambda w^{T} \Sigma w$

みたいなものを解くわけですが、

SBMで解くためには、以下の形（イジング問題）に変換する必要があります。

$x$ を成分が $0$ または $1$ からなるベクトルとして $min_{x} x^{T}Z x$

この $Z$ をSBMのインスタンスに投げることで、答えのresponseが返ってくるイメージです。

問題の変形

さて、上で見た通り通常のポートフォリオ最適化から

SBMで受け取れるように変形する必要があります。

その際、ウェイト $w$ を求める問題から数量 $x$ を求める問題に変わってしまいますが

それは特に問題ないかなと思います。

2進展開

たとえば、12という数字は、以下のように2進展開できます。

$12 = 0 * 2^{0} + 0 * 2^{1} + 1 * 2^{2} + 1 * 2^{3}$

このテクニックを使うことで、

12という数字は $(0,0,1,1)$ という01ベクトルと対応付けることができそうです。

すると、例えば期待リターンが3の証券を12個保有した際の

期待リターンは以下のようにあらわすことができます。

$3 * (1, 2, 4, 8) (0, 0, 1, 1)^{T} = 36$

さらに、複数証券のときは、 $d = (1,2,4,8)$ として

期待リターンが3の証券を12個、3.5の証券を10個保有した際は

$(3, 3.5)C(0,0,1,1,0,1,0,1)^{T}=71$

ただし $C$ は対角成分に $d$ が並んでいる2*8行列

のようにあらわすことができそうです。

このように考えることで、

$\mu^{T} C x$

という形で、証券を $x$ だけ保有した際の期待リターンが求まる。

ただし、これではまだ二次形式になっていないのでSBMには投げれない。

そこで、01ベクトルはその成分を2乗しても変わらないことを利用して

$\mu^{T} C x = x^{T} diag(\mu) \otimes diag(d) x$

とすることで、はれて二次形式となる。

そして、リスクについては

$x^{T}C^{T} \Sigma C x$

とすれば、 $x$ だけ保有した際のリスクになるので

$x^{T} (C^{T} \Sigma C - diag(\mu) \otimes diag(d)) x$

を最小化することで、いわゆるポートフォリオ最適化と類似の問題に帰着する。

コード

最小化するべき二次形式の真ん中の行列をまずは計算する。

サンプルとして以下のようなものを想定する。

import numpy as np


Z = np.zeros((2,2))
Z[0][0] = 10
Z[0][1] = 3
Z[1][0] = 4
Z[1][1] = 7

次にこれを所定の形式で吐き出す。

from scipy.io import mmwrite
from scipy.sparse import csr_matrix


mmwrite('Z.mtx', csr_matrix(Z), field='real', symmetry='symmetric')

最後にこれをSBMのインスタンスに投げつけて、

レスポンスを読み取る

from urllib import request


res = request.urlopen(
    request.Request(url='***', data=open('Z.mtx', 'rb'), headers={'Content-Type': 'application/octet-stream'}
)
result = json.load(res.read())

となります。

とはいえ

まぁポートフォリオ最適化だと

整数計画に精緻化するメリットが希薄なのと、

変数の数が2進展開で爆増してしまって

結局スピード面の優位性がそがれるので、

あまり良い事例ではないように感じられます。

東芝社の論文内で述べられている為替の三角裁定や

ダルマキャピタルとの共同研究の株 stat arbのほうが筋が良さそうですね。

ではまた！

今回の関連図書

リンク

そこそこわかりやすく記載されていた気がします。