MATLAB、Pythonで株価予測【バックナンバー】

MATLAB、Pythonで株価予測【バックナンバー】MATLAB/Simulink
MATLAB、Pythonで株価予測【バックナンバー】
スポンサーリンク

はじめに

MATLAB、Pythonを使って株価予測を使用と考えるシリーズ。
と言っても基本的にはフーリエ変換が中心のネタとなる。

そういう意味では予測と言うよりも分析に近いが、
最終的には予測したいんで分析してるって感じでもある。

尚、株価予測に関してはLSTMによる予測もやっているので興味ある方は以下よりどうぞ。

スポンサーリンク

書籍とか

まぁここら辺の書籍は巷に溢れているので適当に選んで読んでいってもOKだと思います。

スポンサーリンク

導入編

  • FFTで振動解析を行うことが多い。
    • 自動車だと静粛性評価などが代表的。
  • 株価も振動っぽいからFFTで解析できるかも?
    • 周波数成分を見れることは間違いない。
    • ただし、そこから予測に至れるかは別問題。
    • ここでの話を鵜呑みにして株売買をしてもそれは自己責任。
  • 業務でFFTを使っていても、FFTそのものが何か知らない人も多い。
    • 見るべき、比較すべきデータが揃っていると割と知らなくても平気。
      • これ自体は標準化、過去データ利用の結果なので褒められるべき事象。
  • FFTを知るには最低限以下の知識が必要。
    • フーリエ変換、DFT、FFTとそれらの逆変換。
  • フーリエ変換について簡単に説明。
    • 実際には逆フーリエ変換も含めないと全体像は見えない。
  • フーリエ変換自体は実用されてるツールなので、数学者と言うよりエンジニア側の領域
  • 逆フーリエ変換を知らずにフーリエ変換だけの性質を見ると不透明な感じがして恐怖感が芽生えている可能性が高い。
スポンサーリンク

FFTへ至る道

  • 逆フーリエ変換が正しいのかフーリエ逆変換が正しいのか。
    • どっちも正しいと思っておいた方が良さそう。
    • 英単語の並びを重視するか、逆変換という日本語としての意味を重視するか。
  • フーリエ変換/逆変換はバンドパスフィルタ利用が有名。
  • フーリエ変換、逆フーリエ変換にバリエーションがある点に注意。
  • フーリエ変換、逆フーリエ変換の元ネタがフーリエの積分公式。
    • f(t)とf(x)は同じものだが、複素指数関数との畳み込み積分を経由しても等しい状態を作れることを示している。
    • 複素指数関数はオイラーの公式より三角関数に展開可能。
      • 畳み込み積分は三角関数とf(t)の内積を示しており、同一角周波数のみが取り出せる理屈。
  • フーリエの積分公式は「とある関数を畳み込み積分を経ても同じ関数に戻せる」と証明されているもの。
    • 複素フーリエ級数、複素フーリエ係数で証明可能だが、ここでは省略。
  • フーリエの積分公式の一部をフーリエ変換と定義した。
  • フーリエ変換の式をフーリエの積分公式に戻すことで逆フーリエ変換の式が完成。
  • フーリエ変換とフーリエ逆変換のバリエーション自体もバリエーションがある。
    • 角周波数表現と周波数表現によるバリエーション。
    • 数式対称性によるバリエーション。
  • 角周波数表現は前回&今回再掲したもの。
    • 周波数表現は角周波数を単純に周波数の式を代入したもの。
  • 角周波数は角速度のスカラー量。
  • フーリエ変換とフーリエ逆変換のもう一つのバリエーションである、数式対称性について。
    • 3パターンある。
      • 1/2πをどちらが持つかって違い。
      • 1/√2πのように折半するパターンもある。
  • バリエーションを認識していないと異なるバリエーションの変換/逆変換の組み合わせを使用してしまい、元の波形に戻らない事故が発生する。
  • 大雑把に逆離散フーリエ変換と離散フーリエ変換の導出を説明。
    • サンプリング都合で積分範囲をプラス側へ。
    • 逆変換を想定した逆行列都合でサンプリング数と級数の数を合わせて正方行列が作れる状態にしておく。
  • 逆離散フーリエ変換の行列演算形式に対してエルミート転置を利用して逆変換を求める。
    • これが離散フーリエ変換になる。
  • 離散フーリエ変換、逆離散フーリエ変換のバリエーションについて。
    • フーリエ変換、逆フーリエ変換の時と同じく数式対称性によるバリエーション。
    • 1/Nをどちらによせるか、折半するか。
    • よって、対になってる離散フーリエ変換、逆フーリエ変換を使用する必要がある。
  • FFT、IFFTの数式上のバリエーションはDFT、IDFTと一緒。
    • 元にしている数式自体は同一。
    • FFT、IFFTはバタフライ変換による高速化を行ってる点で異なるのみ。
  • バタフライ変換を理解するためには回転因子のイメージが重要。
    • オイラーの公式のおかげで複素指数関数と三角関数が紐づく。
  • タイトル詐欺にならないようにMATLAB、Pythonを使って各回転因子を算出して見た。
    • 1/√2のような計算結果にはならないので1/√2=0.7071を想定した見方になる。
    • 前回の回転因子と同一の結果が得られた。
  • 虚数表現はMATLABはi、Pythonはjとなっている。
  • 回転因子を元に行列表現してみた。
    • 回転位置による最適化が可能。
    • 必ず実数になる点が存在することによる最適化が可能。
    • 対角線による最適化が可能。
  • ここまででもかなり便利ではあるが、さらにバタフライ演算をするための最適化もある。
  • Nの回転因子にN/2の回転因子を含めることが可能
  • 複数段の行列に分解可能。
    • つまり、演算を分解できる。
  • 最終的に残る値はかなり限られる。
    • これを利用してバタフライ演算を行うことになる
  • バタフライ演算を図示した。
    • 演算が交差している様が蝶々のようなのでバタフライ演算と呼ばれている。
  • 入力サンプリングの入れ替えルールはビットリバースに準じている。
  • これら一連の流れを「CooleyTukey型FFTアルゴリズム」と呼ぶ。
    • このアルゴリズムは入力サンプリングが2のべき乗であることが前提。
スポンサーリンク

ETF/VTI

  • やっと株の話に突入
  • 採用する商品はVTIという米国ETF
    • 理由としては、「外乱を可能な限り減らし、分析難易度を下げたい。」
    • それでも外乱は0ではないが、コロナ禍、某戦争、円安、円高など比較的分かり易いものとなる。
      • 個別銘柄だと、その市場や競合他社動向など調べるべきことが細かくなりすぎる。
  • VTIについて少し掘り下げを開始。
  • ベンチマークは「CRSP USトータル・マーケット・インデックス」。
    • 米国株の4000銘柄の変動。
  • バンガード社は他の指標を元にした商品も出している。
    • VTIとよく比較されるのがVOO。
    • VOOはS&P500という米国株500銘柄と連動。
  • シミュレーションを恐らくやると思うので、その前提条件を揃えるため手数料とか税金について調査。
    • SBI証券の情報をベースに調査。
  • 手数料。
    • 買付手数料\0
    • 売却手数料0.45%
  • 税金。
    • 譲渡所得。
      • 合計で約20%
        • 所得税15%。
        • 住民税5%。
    • NISA枠は一旦無視する。
  • VTIチャートを取得する必要がある。
    • 各証券会社だけでなくGoogle検索もでチャート自体はすぐ見れる。
  • 分析する上でcsv形式でチャートが公開されてるとうれしい。
    • Yahooファイナンスの英語サイトでcvs形式で公開されている。
    • この手のURLは変わる可能性が高いので、探し方の一連の操作を確認。
  • チャート取得範囲は2020年10月~2021年9月とする。
    • 直近の米国事業年度。
  • VTIチャート各列を確認。
    • 各列
      • Open(始値)
      • High(最高値)
      • Low(最安値)
      • Close(終値)
      • Adj Close(調整後終値)
      • Volume(売買高)
    • Close(終値)を元に分析、予測を実施する予定。
  • データ整形方針。
    • 単純に平均値を出して、各要素から平均値を引く。
    • バイアス除去になり、周波数0[Hz]が無くなる。
    • MATLABのVersionによってはcvs読み込みに難点があるので、可能な限りシンプルなデータに整形しておく。
  • グラフ表示してみた。
    • ずっと右肩上がり。
    • 本来であれば、一回投資したら放置が正解。
スポンサーリンク

FFT/IFFT

  • FFT、IFFTの入出力って実は良く分かってない。
    • よって、自明且つシンプルな波形を入れて評価して見た方が良い。
      • 自明且つシンプルな波形はsin波とかそれらの合成波。
  • 実験はMATLABで実施するが、Pythonコードを起こす予定。
    • ぶっちゃけメンドクサイとは思ってるけど頑張ってやる。
  • sin波でFFT、IFFTを実施。
    • パッと見ちゃんと元に戻ってるのは確認。
  • FFTの出力である周波数の分布は前半と後半で意味が異なる。
    • 後半が前半の複素共役に当たり、IFFT時に虚数部を相殺する役割を追っている。
    • 特定周波数を取り出す場合は複素共役部分も一緒に取り出す必要がある。
  • 引き続きsin波をFFTに入れる実験継続。
  • 今回は入力期間を2πから4πに増やしてみた。
    • これにより、結果的に期間内のsin(x)の振動は増える。
    • よって、FFTの結果としてのsin(x)の周波数は1Hzではなく2Hzとなる。
    • あくまで、入力サンプリングを1周期とした周波数である点に注意。
      • 物理的な周波数とは異なる。
  • これまでMATLABで実験してきたので、Python版コードも作成。
    • 結果は同一と見なせる。
  • MATLABとPython(Numpy)のFFT、IFFTは同一の数式を元にしている。
    • よって、互換性ありと見なしてOKそう。
    • 演算誤差の方が異なるが無視してもOKなレベル。
  • 複数の周波数のsin波を合成したもの対してFFT&IFFT実施。
    • 波形の合成は単純に足し算するだけ。
    • 想定通りの周波数分布になった。
  • 試しに入力サンプリング期間を2倍に伸ばしてみた
    • 想定通り、周波数が2倍になる分布に変化。
  • 複数の周波数のsin波を合成したもの大してFFT&IFFT実施のPython版。
    • 当然ではあるが、同一の結果が得られた。
  • 加えて、入力サンプリング期間を延ばしたものも実施。
    • これも当然、同一の結果が得られた。
  • FFT,IFFTの理屈は兎も角として使い方に関してはなんとなく慣れてきたところ。
  • IFFTにも活躍してもらうため簡易的なバンドパスフィルタを実施予定。
  • 周波数分布関数の後半に複素共役が居るため、これも同等の処置が必要。
  • 複素共役の位置を分かりやすくするため、マイナス側に持ってくる予定。
    • ベクトルに対するシフト、ローテーションで対応可能。なはず。
  • MATLABでベクトルローテーションをさせたい場合はcircshift関数を使えばOK。
    • ただし、次元指定の罠がある。
      • デフォルト次元が列方向なので、行ベクトルに対して行うとローテーションされない。
      • (実際には列方向にローテーションはされてると思う)
      • [0,3]のように列行それぞれにローテーション数を指定する。
  • FFT出力の周波数分布をローテーションをMATLABで実施。
    • 複素共役が0点を中心とした線対称になるように配置。
    • この配置の方が確認し易さ、処理のし易さが増す想定。
      • よって、配置が意識できていればやらなくてもOK。
      • 効能は実際の処理をする際に確認。
  • Python(Numpy)によるベクトルのローテーションをお試し。
    • 無事ローテーション可能。
  • Pythonはベクトルに種類がある。
    • ベクトル。
      • 転置等の行列由来の演算ができない。
    • 行ベクトル。
    • 列ベクトル。
    • 行列演算を意識する際はベクトル以外の定義をする必要がある。
    • ※ 今回は「ベクトル」でOK。
  • Python(Numpy)でFFT出力の周波数分布をローテーションするコードを作成。
  • 上記の動作確認。
    • MATLABの結果と同じく、負の周波数側に複素共役を持ってこれた。
  • これで超簡易バンドパスフィルタの事前準備が整った。
  • MATLABで超簡易バンドパスフィルタ実施。
    • コード開示&結果確認。
    • 想定通り、3Hzだけ抽出で来た。
  • 3[Hz]と-3[Hz]以外を0にしてる部分のコードが妙。
    • 論理インデックス検索という手法を使っている。
  • MATLABのインデックス検索には線形インデックス検索と論理インデックス検索がある。
  • 線形インデックス検索はC言語の配列の添え字の考え方と一緒だが、添え字に設定できる値がベクトルにできる。
    • C言語はスカラーのみ。
  • 論理インデックス検索は渡すベクトルの1かtrueのところだけが参照できる。
  • それぞれ書き換えも可能。
  • Python(Numpy)による論理インデックス検索と線形インデックス検索が可能か確認。
    • 一応、可能そう。
    • ただし、行列に対する線形インデックス検索はちゃんと行と列を指定する必要あり。
    • MATLABと同様にするにはrehapeでベクトルに直す必要がある。
    • あと、MATLABは1オリジン、Pythonは0オリジンな点も注意。
  • Python(Numpy)によるバンドパスフィルタのコード作成。
  • 上記コードを実行して見た。
    • MTALABと同じ結果が得られることを確認。
    • よって、MATLAB、Python双方でFFT、IFFTによる特性周波数の抽出が可能と言える。
スポンサーリンク

周波数解析

  • FFT,IFFTの使い方及び、特定周波数の抽出方法が分かったところで今後の方針を考える。
    • とりあえずVTIチャートにFFTかけて周波数特性見て見る。
    • 後のことは見てから考える。
  • 大雑把すぎる方針だが、ホントやってみないとなんもわからん。
  • FFTの出力をローテーションしている都合、サンプリング数は偶数が望ましい。
    • 頑張って調整すれば奇数でも行けるはずだがメンドイのでやらない。
  • 無事、VTIチャートの周波数特性及びそこからのVTIチャートへの逆変換ができた。
    • しかし、VTIチャートの周波数特性自体に何かしら問題が・・・。
  • Python(Numpy)でVTIチャートにFFT、IFFTをブチかます。
    • MATLABと同一の結果が得られた。
    • これによりMATLABとPython(Numpy)の両方で同じレベルで実験が進められる。
  • Numpyにcsv読み込みの機能がある。
    • 他にも手段はあるが、今回はNumpyで実施。
  • 今回使用するVTIチャートの問題について考察。
    • あまりにもシンプル過ぎて、分析し易い特徴が捕まえられない可能性あり。
    • 一応5[Hz]が若干飛び出ているが、特徴になり得るかは疑問。
    • とりえあえず、現状のチャートでやってみて、分析が難しいようであれば再度方針検討。
  • VTIチャートの周波数特性から5[Hz]と抽出してみる方針に。
  • 抽出した5[Hz]は元のVTIチャートの振幅と比べると遥かに小さい可能性が高い。
    • よって、グラフで比較する際はIFFT側の結果を増幅してあげた方が良い。
  • とりあえず、現状のVTIチャートに対してMATLABで5Hzを抽出するコードを作成。
    • IFFT後の波形の振幅の増幅は最大値を比較して、その比率を使用して増幅係数を決定。
  • 一応5Hzの抽出はできたが・・・。
    • ここらへんの考察はPythonコードを作成した後に実施予定。
  • 前回MATLABで作ったVTIチャートから5Hzを抽出するコードのPython版を作成。
    • 振幅調整も同じ処理で対応。
  • MATLABと同じ結果が得られたことは確認。
  • 問題は、これから何を分析できるかと言う点だが、そこは次回。
  • VTIチャートとIFFTの結果に対して考察。
    • 5Hzではさほど何かを示しているデータには見えない。
  • 試しに13Hzを中心として12Hz~14Hzを抽出。
    • こちらは細かい山と谷を捕まえて居そう。
    • 期間が短いことから値動きも小さい。
  • 今後の方針としては2022年も含めたVTIチャートで試してみる。
  • 新しく2021年6月から2022年5月のVTIチャートを取得。
    • MATLABでplotしてみた。
    • いい感じに乱高下している。
  • MATLABとPython(Numpy)のFFT、IFFTで元の波形に戻せるかを確認。
    • 共に戻せることが確認できたので、周波数解析できそう。
  • 新VTIチャートの周波数特性を確認。
    • 1[Hz]、3[Hz]、5[Hz]、7[Hz]あたりが突出している。
  • まずはお試しで3[Hz]を抽出して確認。
    • かなり特徴を表していることが分かる。
    • Python(Numpy)でも同様のことができることを確認。
  • 新VTIチャートから各種周波数を抽出して比較。
    • 3[Hz]が筋が良さそう。
    • ついで5[Hz]。
    • 7[Hz]は細かい特性は掴んでいるが、ちょっとイマイチ。
  • 今回の情報だけを見ると3[Hz]で売り買いすればOKということになるが・・・。
  • 新VTIチャートからの抽出周波数を1個ではなく、複数にすることで筋の良さそうな特性が出てくる。
  • しかし増やせば良いというものではない。
    • 増やせば増やすほど元の新VTIチャートに近づくだけ。
      • 単にローパスフィルタを掛けただけになることも。
    • 何パターンか出してみて自分自身が信用できそうなものを探すって流れになる。

スポンサーリンク

収支シミュレーション

  • 収支シミュレーションをする上で、今後の方針を決めた。
  • 基本的にはプログラム的に算出するが、実際の売買時の計算は手動で。
    • 気が向いたら自動演算化するかも。
  • 極大値、極小値を特定する必要があるが、これはそれほど難しくはない。(つもり)
  • いままではExcelで終値から平均値を引いた値と算出したものをcsvにしていたが、MATLAB、Pythonで終値を取り込んでから平均値を引く方式に変更。
  • MATLAB、Pythonともにmeanという関数/メソッドで平均値算出可能。
    • Pythonはaverageという加重平均を算出するメソッドが存在。
  • 微分して0になれば極値。
    • 極大値は、「微分値がプラス→0→マイナス」となるところ。
    • 極小値は、「微分値がマイナス→0→プラス」となるところ。
  • 「微分して0」でも極値にならないパターンもある。
    • 3次関数とかが代表的。
  • 極大値と極小値の特定のMATLABコードを作成。
  • 上記コードの動作確認。
    • 極大値に赤丸、極小値に青丸を置いてる。
  • 一部問題点あり。
    • 最初に極小値が来ることを想定している。
    • しかし、最初に買付をする想定なので、むしろ今回のコードの方が都合が良い。
  • 極大値と極小値の特定するコードのPython(Numpy)版を作成。
  • 動作としては反転波形も含めてMATLABと同一。
  • 実際にはこれらコードにく分けて以下が必要だが、本番コード作成時に盛り込む。
    • 各プロット時の値の取得。
    • ドルから円へ変換。
      • $1=\127で計算する予定。
  • 売却、買付タイミング時のVTI単価特定コード(MATLAB版)を作成。
  • ついでにバンドパスフィルタの部分をちょい改修。
  • Python(Numpy)側のコードも似たような感じで改修予定。
  • 売却、買付タイミング時のVTI単価特定コード(MATLAB版)の実行結果を確認。
    • グラフで確認。
    • 拡大グラフで確認。
    • コマンドウィンドウ出力を確認。
      • 一部を除いて利益が出そうな数値にはなっている。
  • 売却、買付タイミング時のVTI単価特定コードのPython版を作成。
  • MATLABの時に実施した論理インデックス検索のやり方を修正。
    • このために複素共役を負の周波数側に持ってきた。(忘れてたけど)
  • Python版の売却、買付タイミング時のVTI単価特定コードを実行。
    • グラフはOK。
    • VTI単価出力もOK。
  • 極大値、極小値のインデックスがMATLABで実施した時と異なるが、これはオリジンのせい。
    • MATLABは1オリジン、Pythonは0オリジン。
  • 収支シミュレーションの方針を確認。
    • 買付手数料:買付金額の0%。
    • 売却手数料:売却金額の0.45%。
    • 税金:利益の20%。
    • 売買単位:1口。
    • 買付余力:100万円。
  • VTIなど海外ETFは1口から売買可能。
    • というか100口、100株が売買単位なのは日本特有の文化。
  • ついに収支シミュレーション実施。
    • 収益としては11万円弱。
  • イマイチな結果にも見えるが、最大単価時の一回の売買よりも収益は出ている状態。
    • たまたまの可能性は高い。
  • 前回の収支シミュレーションの考察を実施。
    • 明らかに損がでる売買も行っている。
    • しかし、そこで売買しないと買付余力がなくなり、機会損失になる場合もある。
    • 買付余力があれば、損が出る売買は見送れるので買付余力を如何に残せているかも重要にはなってくる。
  • 同じシミュレーションを個別株で行えないか?
    • 基本ロジックは一緒なので可能。

スポンサーリンク

個別株シミュレーション

  • 分析予測対象の個別株を選んできた。
    • 大企業株のため、VTIの特性に似てはいるが、細かい上下が見て取れる。
  • とりあえず周波数解析実施。
    • 10[Hz]あたりが突出している。
  • 9~11[Hz]の範囲だけを残してIFFT実施。
    • 筋が良さそうなのでこれをベースに掘り下げる予定。
  • 10[Hz]に於ける売買タイミングと株価を確認。
    • 10回の売買が発生する。
  • 収支シミュレーションを行う上でのパラメータ決め。
    • VTIの時の差分は以下。
    • 売買手数料は0円。
    • 売買単位は10株
    • 買付余力は75万。
      • 100万以下の売買を前提として売買手数料を0円にしているため。

スポンサーリンク

コメント

  1. あーる より:

    お聞きしたいことがあります。
    フーリエ変換について、①∫f(x)e^(-2πxξ)dxについてですが、eの指数の中のxのみx+C(0≦C<1/ξ)に置き換えて(②)、Cを変化させながらフーリエ変換を計算した場合。出力はCの値に応じてある程度変化が得られるのでしょうか?

    また、①式と②式、また、全てのCの値についての②式の出力の絶対値を全て積分した出力(③)を実際の音波に作用させて、出力を周波数毎の振幅を表すグラフとして比較した場合、どれが、元の音波の実際の周波数振幅スペクトルに近いと言えるのでしょうか?

    • KEIKEI より:

      ①②に関してですが、
      「時間シフトの法則」が該当するのではないでしょうか?
      よって、
      {∫f(x)e^(-2πixξ)dx}e^(-2πiC)

      ③については、ちょっと状況が見えておらず、回答が難しいです。

      • あーる より:

        ありがとうございます。今、人間の声の周波数成分を簡単に分析したくて、Pythonでフーリエ変換を使い、周波数毎の振幅を計算するためのプログラムを作っていました。

        ③については、サンプリングの開始時間を細かくずらした上で、サンプリング、フーリエ変換して、和を取ったら、sin(2πxξ)の0から2πの位相に対して均等に波を拾うことができるかな〜などと考えて、質問しました。

        質問をするために自分でも平行移動の法則を事前に調べてチェックしていたのですが、もし平行移動により出力の数値が変動するのでしたら、サンプリングの開始時間をずらすと微妙に出力にゆらぎが出るのかもしれない…と思いお聞きできる人を探していました。いきなりお聞きして申し訳ないのですが、もし良ければご意見教えて頂けると嬉しいです。

        • あーる より:

          そっか、回答は難しいと書いてありますね。ここまでで大丈夫です。お邪魔しました。回答頂きありがとうございました。

          • あーる より:

            簡素な返信になってしまい申し訳ありません。。

            感謝しております。多謝。

          • KEIKEI より:

            お力になれず申し訳ないです。

            ③の答えというわけではないですが、
            やりたいこととしては「ケプストラム分析」が近いかもしれません。
            周波数領域に変換後にローパスフィルタをかけるイメージで、周波数特性の平滑にはなります。

タイトルとURLをコピーしました