JP2019125306A - データ処理方法、データ処理装置およびプログラム - Google Patents

Abstract

【課題】高頻度非同期データの分散を求める。【解決手段】データ処理方法のステップ１０は時系列データを取得する。ステップ２０は、時系列データの時間変化の要因である説明変数を取得する。ステップ３０は、時系列データと説明変数とに基づいて、疑似対数尤度関数を最大化することにより、時系列データの分散を決めるパラメータを最適化する。ステップ４０は、最適化されたパラメータに基づいて、時系列データの分散を推定する。【選択図】図２

Description

本発明は、データ処理方法、データ処理装置およびプログラムに関する。

高頻度かつ非同期に発生する時系列データを統計的に扱う手法として、パラメトリック・モデル下における最尤型推定量に関する研究がされている（例えば、非特許文献１参照）。また時系列データを取り扱うニューラル・ネットワークとして、リカーレント・ニューラル・ネットワーク（Ｒｅｃｕｒｒｅｎｔ Ｎｅｕｒａｌ Ｎｅｔｗｏｒｋ。以下、「ＲＮＮ」という）の研究がされている。

"Ｐａｒａｍｅｔｒｉｃ Ｉｎｆｅｒｅｎｃｅ ｆｏｒ Ｎｏｎｓｙｎｃｈｒｏｎｏｕｓｌｙ Ｏｂｓｅｒｖｅｄ Ｄｉｆｆｕｓｉｏｎ Ｐｒｏｃｅｓｓｅｓ ｉｎ ｔｈｅ Ｐｒｅｓｅｎｃｅ ｏｆ Ｍａｒｋｅｔ Ｍｉｃｒｏｓｔｒｕｃｔｕｒｅ Ｎｏｉｓｅ"、Ｏｇｉｈａｒａ Ｔｅｐｐｅｉ、Ｔｈｅ Ｉｎｓｔｉｔｕｔｅ ｏｆ Ｓｔａｔｉｓｔｉｃａｌ Ｍａｔｈｅｍａｔｉｃｓ、Ｂｅｒｎｏｕｌｌｉ、ｉｎ ｐｒｅｓｓ．ｈｔｔｐｓ：／／ｗｗｗ．ｅ−ｐｕｂｌｉｃａｔｉｏｎｓ．ｏｒｇ／ｉｍｓ／ｓｕｂｍｉｓｓｉｏｎ／ＢＥＪ／ｕｓｅｒ／ｓｕｂｍｉｓｓｉｏｎＦｉｌｅ／２４９４８？ｃｏｎｆｉｒｍ＝ｆ４５ｃｅ６４０

観測された時系列データを基にその分散や共分散を計算することは、様々な局面で有用なデータ処理である。ここで分散とは、一般にデータの散らばり度合いを表す値であり、観測されたデータと平均値との差（偏差）の二乗の平均として求められる。また共分散とは、一般に２種類の観測データ同士の関係を表す値であり、各観測データに関する偏差同士の積の平均として求められる。

例えば株式資産のリスク管理において、株価の時系列データの分散・共分散を特定することは極めて重要である。すなわち分散は、当該銘柄の株価変動リスク（ボラティリティ）や値動きの荒さを反映した数値であり、これを把握することはポートフォリオの運用に欠かせない。また共分散は、異なる銘柄同士の株価変動の連動性を反映した数値であるため、やはりその傾向を把握することは重要である。従来、株価の分散・共分散の把握は、その日の終値など１日１回乃至数回程度観測される低頻度のデータを用いて行われてきた。これに対して、時々刻々変動する日内のすべての取引に関する株価データは、株価変動に関する多くの情報を含むため、証券市場の構造分析などに非常に有効なものであると考えられる。

ここで、日内のすべての取引で観測される株価のようなデータには、他の一般的な時系列データ（特に低頻度データ）とは異なる際立った特徴がある点に留意する必要がある。１つはその高頻度性である。株価データは、当該銘柄の証券の取引が成立した瞬間に発生するものである。表１は、ある自動車会社関連株に関し、約３０秒の間に実際に発生した証券取引の状況を表す。表１に示されるようにその発生頻度は通常、数秒または数百ミリ秒に１回である。
これは１日の取引時間全体（例えば５時間）に対して、極めて高い頻度であるといえる。今１つの特徴は観測データ同士の非同期性である。一般に異なる銘柄の証券が取引されるタイミングは異なるため、これらの株価データの発生時刻は一致しない。すなわち、異なる銘柄の株価データの観測は非同期的になされる。

株価のような高頻度かつ非同期で発生するデータ（以下、「高頻度非同期データ」と呼ぶ）には、データの膨大さに加えて、分散・共分散計算などの統計解析を困難にする２つの特徴的な要因が伴うことが分かってきた。１つは「マイクロストラクチャー・ノイズ」と呼ばれる仮想的な観測誤差であり、今１つは「非同期観測」に伴う問題である。以下、この２つの要因について説明する。

（マイクロストラクチャー・ノイズ）
株価のボラティリティの推定値と考えられる値の１つに、実現ボラティリティ（Ｒｅａｌｉｚｅｄ Ｖｏｌａｔｉｌｉｔｙ。以下「ＲＶ」と呼ぶ）がある。ＲＶは、時間的に隣接する株価（または時間的に隣接する株価の対数値）間の差の２乗値を、一日の取引時間内で合計して得られる値であり、観測データを｛Ａ_ｔｋ｝^Ｌ _ｋ＝０とすると以下の式で定義される。
株価の対数値（以下、「対数株価」という）Ｙ_ｔが後述の拡散過程に従い、｛Ｙ_ｔｋ｝^Ｌ _ｋ＝０が観測されると仮定した場合、ＲＶは、ｍａｘ_ｋ（ｔ_ｋ−ｔ_ｋ−１）→０の極限（高頻度観測極限）でＹのボラティリティ＜Ｙ＞_１に収束する。

しかしながら実証研究では、高頻度観測極限でもＲＶが収束せず、逆にデータが高頻度になるとともにＲＶが急激に増加する現象が確認されている。図１は、表１に記載された株価の実データを基に、隣接する時間間隔を３０秒から１８００秒まで変えて計算したＲＶのプロットである。図１に示される通り、ＲＶは、時間間隔（頻度）が３００秒以上では０．０００１付近で安定しているが、時間間隔が３００秒以下になると急激に増加して発散する。このような、観測頻度を高くするとＲＶが急激に増加する現象、すなわち観測頻度の増加とともにＲＶが過大に評価される現象は、実際の株価データの観測時に仮想的な観測ノイズが混入することが原因であると解釈されている。高頻度観測データを解析する際に現れるこの特有のノイズは、「マイクロストラクチャー・ノイズ」と呼ばれている。

（非同期観測）
前述のように、ある銘柄の株価データが観測されるのは、当該銘柄の証券取引が成立した瞬間である。従って複数銘柄の株価データの観測時刻は一致しない。すなわち異なる銘柄の株価データは「非同期に観測される」。このような２種類の非同期に観測される時系列データに対して、線形補間や直前データによる補完などといった簡易な同期化を用いて共分散を算出すると、共分散の深刻な過小評価を招くという問題が発生する。

以上述べたように、高頻度非同期データでは、マイクロストラクチャー・ノイズや非同期観測といった観測上の複雑さの存在が、その統計解析を困難なものとしている。このような観測上の複雑に対処するための先行技術として、パラメトリック・モデル下における最尤型推定量を利用するものがある。これは、対数株価の時間変化が未知パラメータを含む確率微分方程式を満たすと仮定して、疑似尤度関数を用いてこの未知パラメータを推定することによりＲＶを算出するものである（非特許文献１参照）。この手法は、高頻度非同期データに伴う観測上の複雑さに一定程度対応できるものである。

しかしながら株価データの場合、マイクロストラクチャー・ノイズや非同期観測といった観測上の複雑さに加えて、モデル自体の複雑さも存在する。これには主に「日中季節性」と「長期記憶性」がある。日中季節性とは、一日の取引開始時・取引終了時に取引が集中するという性質である。また長期記憶性とは、一度ボラティリティが上昇すると、それが長期にわたって継続されやすいという性質である。前述の手法を用いても、このようなモデル自体が複雑なものを包括的に扱うことのできるパラメトリック・モデルを事前に設定することは難しい。

こうした場合の有力なアプローチとして、機械学習を用いた手法が考えられる。すなわち、事前にモデルを設定することに代えて、観測されたデータから当該データのモデル構造を学習させるというものである。機械学習は、大量の学習データから当該データの特性を学習させ、コンピュータを用いた計算から背後にあるモデル構造を特定する。近年、時系列データを取り扱うニューラル・ネットワークとして、ＲＮＮなどが研究されている。しかしながら、多くの機械学習はユークリッド空間上のデータに対する手法であるため、観測上の複雑さを持つ高頻度非同期データに直接適用することはできない。

これに対し本発明者は、高頻度非同期データを時間の関数とみて関数空間上の要素として扱い、機械学習理論をこの空間上の理論へと発展させることにより、当該データを確率過程論と融合させて扱うことができることに気が付いた。具体的には、高頻度非同期データを扱うことのできる疑似対数尤度関数をニューラル・ネットワークと融合することで、効率的な学習が実現できると考えられる。しかしながら、疑似対数尤度関数をニューラル・ネットワークにそのまま適用すると、疑似対数尤度関数にバイアスが発生する。このため計算が一定の速さで収束しないという問題がある。これは計算精度の低下を招く。従ってこのバイアスをいかに補正するかが課題となる。さらに疑似対数尤度関数の計算では、逆行列を算出するときの計算負荷が非常に大きくなる。従って、この逆行列計算を高速化することも実用上の大きな課題となる。

上記課題を解決するために、本発明のある態様のデータ処理方法は、時系列データを取得するステップと、時系列データの時間変化の要因である説明変数を取得するステップと、時系列データと説明変数とに基づいて、疑似対数尤度関数を最大化することにより、時系列データの分散を決めるパラメータを最適化するステップと、最適化されたパラメータに基づいて、時系列データの分散を推定するステップと、を備える。

パラメータを最適化するステップは、バイアスを補正するステップをさらに含んでもよい。

パラメータを最適化するステップは、高速近似計算を用いて逆行列を計算するステップをさらに含んでもよい。

時系列データは複数種類の時系列データであってもよく、時系列データの分散を推定するステップは、時系列データの共分散を推定するステップをさらに含んでもよい。

パラメータを最適化するステップは、ニューラル・ネットワークを用いてもよい。

本方法は、説明変数を平均化するステップをさらに備えてもよい。

疑似対数尤度関数を最大化するステップは、高速近似計算を用いて、行列式の対数を計算するステップをさらに含んでもよい。

時系列データは、株価データであってよい。

本発明の別の態様のデータ処理装置は、時系列データを取得する時系列データ取得部と、時系列データの時間変化の要因である説明変数を取得する説明変数取得部と、時系列データと説明変数とに基づいて、疑似対数尤度関数を最大化することにより、時系列データの分散を決めるパラメータを最適化するパラメータ最適化部と、最適化されたパラメータに基づいて、時系列データの分散を推定する分散推定部と、を備える。

本発明のさらに別の態様のプログラムは、時系列データを取得するステップと、時系列データの時間変化の要因である説明変数を取得するステップと、時系列データと説明変数とに基づいて、疑似対数尤度関数を最大化することにより、時系列データの分散を決めるパラメータを最適化するステップと、最適化されたパラメータに基づいて、時系列データの分散を推定するステップと、をコンピュータに実行させる。

なお、以上の構成要素の任意の組合せ、本発明の表現を装置、方法、システム、記録媒体、コンピュータプログラムなどの間で変換したものもまた、本発明の態様として有効である。

本発明によれば、高頻度非同期データの分散・共分散を求めることができる。

実際に観測された株価データに関し、隣接する時間間隔を変化させて計算したときのＲＶ値を示す図である。 一実施形態に係るデータ処理方法の処理フロー図である。 一実施形態に係るデータ処理方法に用いるニューラル・ネットワークを示す図である。 実際の午前の取引で発生した日内株価データを基に予測した株価のボラティリティを示す図である。 実際の午後の取引で発生した日内株価データを基に予測した株価のボラティリティを示す図である。 電力会社株のＲＶ値について、実データに基づく算出値と本手法を用いた予測値との比較を示す図である。

具体的な実施例を説明する前に、高頻度非同期データの例として株価データを取り上げて、基礎となる事項を説明する。以下、記号を簡略化するため時間単位を変換し、観測の最終時刻を１に規格化したデータを取り扱う。
複数銘柄の株価を考えたとき、時刻ｔにおけるｉ番目の銘柄の株価Ｙ^ｉ _ｔは、以下のような確率過程に従って時間変化するものと仮定する。この確率過程は拡散過程と呼ばれる。
ここで、Ｘ_ｔは説明変数、ａ^ｉおよびｂ_ｉｊは未知の非線形関数、Ｗ^ｊ _ｓはブラウン運動である。説明変数は上記の銘柄の株価の時間的変化の要因となる変数であって、例えば当該銘柄の値動き、別銘柄の値動き、日経平均株価等の株価指数、取引量などに関する情報である。上記の拡散過程の式の右辺第１項は時間に比例して動く部分を表し、右辺第２項はランダムな拡散項を表す。高頻度非同期データの場合、この拡散項を求めることが課題となる。右辺第２項に含まれる非線形関数ｂ_ｉｊ（ｔ、Ｘ_ｔ）は、時刻ｔにおける瞬間的なボラティリティ（以下、「瞬間的ボラティリティ」という）を表す量である。ボラティリティ＜Ｙ^ｉ＞_１および共変動＜Ｙ^ｉ、Ｙ^ｊ＞_１は、瞬間的ボラティリティｂ_ｉｊ（ｔ、Ｘ_ｔ）に基づいて以下のように算出される。

（実施例）
以下、本発明に係るデータ処理方法の一実施形態を詳細に説明する。
図２は、本実施形態の処理フロー図である。ステップ１０で本方法は、観測時刻ｓ^１ _１、．．．ｓ^１ _ｍおよびｓ^２ _１、．．．、ｓ^２ _ｎにおけるそれぞれの時系列データＹ^１ _１、．．．Ｙ^１ _ｍおよびＹ^２ _１、．．．、Ｙ^２ _ｎを取得する。ここで、上付き添字はデータの種類を表し、下付き添字はデータを取得した順番を示す。例えば時系列データが銘柄１および２の日内対数株価だった場合、ｓ^１ _ｉは銘柄１のその日のｉ番目の取引発生時刻を表し、Ｙ^１ _ｉは時刻ｓ^１ _ｉにおける銘柄１の対数株価を表す。なお、本例では時系列データの種類が２つの場合について説明するが、これは１種類のみであってもよく、３種類以上であってもよい。時系列データが１種類のみの場合は、後述の分散のみが計算され、共分散は計算されない。

ステップ２０で本方法は、上記の時系列データの説明変数（ｔ、Ｘ_ｔ）を取得する。

ステップ３０で本方法は、時系列データの分散を決めるパラメータを最適化する。以下、瞬間的ボラティリティｂ（ｔ、Ｘ_ｔ）の推定値ｂ^ｌｉ _ｋはパラメータβを含むとして、ｂ^ｌｉ _ｋ（β）で表す。βを最適化することにより、ｂ^ｌｉ _ｋ（β）の最適値が定まり、これによりデータの分散・共分散が計算される。以下に説明するように、βの最適化は、疑似対数尤度関数を最大化することにより達成される。以下、ニューラル・ネットワークを用いてβの最適化を実行する例を説明する。図３は、多次元パラメータβを持つニューラル・ネットワークｂ（χ；β）を示す。この場合、ニューラル・ネットワークｂ（χ；β）への入力値は説明変数（ｔ、Ｘ_ｔ）であり、ニューラル・ネットワークｂ（χ；β）からの出力値は瞬間的ボラティリティの推定値ｂ^ｌｉ _ｋ（β）である。パラメータβ＝｛β^ｋ _ｉ，ｊ｝に対し、中間層の値ｕ^ｊ _ｋを以下のように定める。
（ｕ^０ _１、ｕ^０ _２、ｕ^０ _３）＝（ｔ、Ｘ_１、Ｘ_２）
ｕ^ｊ _ｋ＝（Σ_ｉβ^ｋ _ｉ，ｊｕ^ｋ−１ _ｉ）_＋
ｂ_ｉｊ＝Σ_ｌβ^Ｋ＋１ _{ｌ，２ｉ＋ｊ−２}ｕ^Ｋ _ｌ
ここで記号（・・・）_＋は、括弧内の値と０とを比較し、いずれか大きい方の値を取ることを意味する。

以下、疑似対数尤度関数Ｈ（β）を最大化することによりβを最適化する処理について説明する。
先ず、観測区間全体を等間隔のブロック［ｕ_０、ｕ_１）、［ｕ_１、ｕ_２）、．．．、［ｕ_ｐ−１、ｕ_ｐ）に分割する。例えば時系列データが日内株価データだった場合、観測区間全体は１日の取引時間全体である。
次に観測の増分Δ_ｋを以下の式で表す。
ここでＴはベクトルや行列の転置を表す。ノイズが加算的で平均０の正規分布に従うとみなすことにより、Δ_ｋは元の確率過程と同様に、局所的に多変量正規分布で近似することができる。
次に分散共分散行列Ｓ_ｋ（β、ν_１、ν_２）を以下の式で定義する。
分散共分散行列Ｓ_ｋ（β、ν_１、ν_２）の各成分は以下の通りである。
ｂ^ｌ _ｋ（β）＝（ｂ^ｌｉ _ｋ（β））_ｉ
（Ｍ_ｌ、ｋ）_ｉｊ＝２δ_ｉｊ−δ_{ｉ、ｊ−１}−δ_{ｉ−１、ｊ}
（Ｌ_ｌ、ｋ）_ｉｊ＝δ_ｉｊ（ｓ^ｌ _ｉ−ｓ^ｌ _ｉ−１）
（Ｌ_{１、２、ｋ}）_ｉｊ＝（ｓ^１ _ｉ∧ｓ^２ _ｊ−ｓ^１ _ｉ−１∨ｓ^２ _ｊ−１）_＋
ここでδ_ｉｊはクロネッカーのデルタであり、Ｘ∧ＹはＸとＹの小さい方の値を取り、Ｘ∨ＹはＸとＹの大きい方の値を取ることを意味する。
そして、ノイズ分散の推定量
に対して、疑似対数尤度関数Ｈ（β）を以下のように定義する。

非特許文献１に記載された技術では、パラメータβの最適値として、Ｈ（β）を最大にするときのβを算出する。しかしながら、以下に説明するように、疑似対数尤度関数をニューラル・ネットワークに適用した場合は、疑似対数尤度関数にバイアスが発生する。これは推定値精度の低下を招く。この問題を解決するために、本発明では発生したバイアスを補正する。以下、バイアス補正処理の一例を説明する。

（バイアスの補正）
ステップ３２で本方法は、疑似対数尤度関数のバイアスを補正する。
非特許文献１に記載された技術では、Ｈ（β）を最大にするときのβの値をβ^＃とすると、
「あるβ_＊が一意に存在してｂ（Ｘ_ｔ）＝ｂ（Ｘ_ｔ；β_＊）」 （１）
という条件が成立するときに、ある確率変数Ａに対し、
β^＃−β_＊〜（ｍ＋ｎ）^−１／４Ａ
が成立している。ここで、ｍおよびｎはそれぞれＹ^１ _ｉおよびＹ^２ _ｉの観測数である。β^＃はパラメータ推定量、β_＊は真のパラメータ値である。すなわちこの場合、観測数が増加するとともに、β^＃は観測数の−１／４乗の速さで真の値β_＊に近づいていく。さらにこの収束の速さは、推定量の種類に関わらず、最も速いものであることも証明されていた。

しかしながら上記の条件（１）は、疑似対数尤度関数をニューラル・ネットワークに適用した場合は通常成立しない。これは疑似対数尤度関数にバイアスが発生するためと考えられる。このようにβ^＃が（ｍ＋ｎ）^−１／４の速さで収束しないことは、βの推定精度が低下することを意味する。

本発明者は、以下の手法を用いてバイアスを補正し、補正後の新たな疑似対数尤度関数Ｈ^＆（β）を適用することにより、Ｈ^＆（β）を最大にするときのβの値β^＆が（ｍ＋ｎ）^−１／４の速さで収束することに気が付いた。

先ず関数ｇ（ｘ）およびｇ^ｊ _ｌ、ｋを以下のように定義する。
ｇ（ｘ）＝ｘ∧（１−ｘ）
ｇ^ｊ _ｌ、ｋ＝ｇ（ｌ／ｍ^ｊ _ｋ）
ここでｍ^ｊ _ｋは、データ列を複数のブロックに分割したときの、ｋ番目のブロックにおけるｊ番目の変数の観測数を表す。
さらに、（ｉ、ｊ）成分が以下のように定義される２×２行列Ｂ_ｋを導入する。

さらに以下の関数を導入する。

ここで、バイアス補正された疑似対数尤度関数Ｈ^＆（β）を以下のように定義する。
Ｈ^＆（β）＝Ｈ（β）＋１／２・Σ^ｐ _ｋ＝２（Ｇ_１、ｋ（β）＋Ｇ_２、ｋ（β）＋Ｇ_３，ｋ（β））
本式に示されるように、従来のＨ（β）に補正項（第２項以降）を付加したＨ^＆（β）を導入すると、Ｈ^＆（β）を最大にするときのβの値β^＆は、（ｍ＋ｎ）^−１／４の速さで収束する。従ってこのようにバイアスを補正することにより、疑似対数尤度関数をニューラル・ネットワークに適用した場合も、パラメータβの最適値β^＆を得ることができる。

（逆行列の高速近似計算）
ステップ３４で本方法は、高速近似計算を用いて逆行列を計算する。
Ｈ^＆（β）を計算するためには、分散共分散行列Ｓ_ｋの逆行列
を求める必要がある。通常ニューラル・ネットワークの計算には、ＧＰＵ（Ｇｒａｐｈｉｃｓ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）による並列化計算が用いられる。しかしながら逆行列計算は並列化が難しいため、実用面での大きなネックとなる。本発明者は、前述の分散共分散行列が特殊な形を持つことに着目し、この特殊な形を利用することにより逆行列を高速に近似計算できることに気が付いた。以下、高速近似計算の一例を説明する。

分散共分散行列の逆行列は、
と置くことにより、以下のように行列の無限級数として表現することができる。
この級数は収束が速いため、有限のｒ（例えばｒ＝１０）までの部分和による近似計算で高い精度の計算結果が得られる。この部分和による近似は行列の積に帰着することから、結局はＥ^ｊ _ｋの高速化だけが課題となる。すなわちＥ^ｊ _ｋの定義式から分かるように、Ｅ^ｊ _ｋの算出には逆行列計算が必要であるため、これを高速化する必要がある。以下、このための高速近似計算の一例を説明する。

まず、
とおく。ｎ^ｊ _ｋはその定義から分かるように、複数のブロックを左から順に並べたとき、ｋ番目のブロックにおける一番右側にあるデータＹ^ｊ _ｉの添字ｉの値である。
これらを基に数列ｃ^ｊ、ｋ _ｉを以下のように帰納的に定義する。
（ｉ≧２）
これを用いることにより、Ｅ^ｊ _ｋは以下のように表現することができる。
ｃ^ｊ、ｋ _ｉは逐次的に計算できることから、結局Ｅ^ｊ _ｋの逆行列計算は、計算負荷の小さい行列の積に計算に帰着する。このようにして逆行列計算を高速化することができる。

（行列式の対数の高速近似計算）
Ｈ^＆（β）の計算、すなわちＨ（β）の計算では、ｌｏｇ ｄｅｔＳ_ｋを算出する必要がある。この行列式の対数の計算も、計算負荷の面で課題となる。本発明者は、以下の手法を用いることにより、これも高速化することができることに気が付いた。
ｌｏｇ ｄｅｔＳ_ｋは、以下のように表現することができる。
これと前述のＥ^ｊ _ｋの表現
から、
となることを用いれば、ｌｏｇ ｄｅｔＳ_ｋを計算することができる。
無限級数
も前述と同様に収束が速いため、有限のｒまでの部分和による近似計算で高い精度の計算結果が得られる。
このようにして行列式の対数の計算を高速化することができる。

パラメータβの最適化は、−Ｈ^＆（β）をニューラル・ネットワークの損失関数としてこれを最小化することによって、換言すればＨ^＆（β）を最大化することによって実現することができる。具体的な計算には、既存のニューラル・ネットワークの最適化手法が適用されてよい。一例として、ＡｄａＤｅｌｔａによる手法が用いられてよい。

ステップ４０で本方法は、ステップ３０で最適化されたパラメータに基づき、時系列データの分散・共分散を計算する。前述のように、パラメータβの最適値が定まると、
ｂ_ｉｊ（β）＝Σ_ｌβ^Ｋ＋１ _{ｌ，２ｉ＋ｊ−２}ｕ^Ｋ _ｌ
から瞬間的ボラティリティの最適値ｂ_ｉｊ（β）が定まる。
これに基づき、Ｙ^ｉのボラティリティ＜Ｙ^ｉ＞_１と、Ｙ^１、Ｙ^２の共変動＜Ｙ^１、Ｙ^２＞_１を以下のように推定することができる。
以上で本実施形態の説明を終える。
なお本発明は高頻度非同期データを対象とする場合に特に有用であるが、低頻度に発生するデータや同期的に観測されるデータの統計解析に対しても十分な効果を奏することはいうまでもない。

（説明変数の平均化）
一般に時系列データが高頻度非同期データだった場合、その説明変数の観測値Ｘ_ｔ＝（Ｘ^ｊ _ｔ）^ｑ _ｊ＝１にもマイクロストラクチャー・ノイズが混入している可能性がある。説明変数に含まれるマイクロストラクチャー・ノイズを除去するために、選択的に、説明変数の平均化処理を行ってもよい。以下、説明変数の平均化処理の一例を説明する。

平均化処理では、観測区間全体を等間隔のブロック［ｕ_０、ｕ_１）、［ｕ_１、ｕ_２）、．．．、［ｕ_ｐ−１、ｕ_ｐ）に分割した上で、ブロック［ｕ_ｋ−１、ｕ_ｋ）における説明変数のＸ^ｉ _ｊの平均値Ｚ^ｊ _ｋを
で定義する。
ここで、ｈ^ｊ _ｋはブロック［ｕ_ｋ−１、ｕ_ｋ）内のＸ^ｊ _ｉの個数、ｔ^ｊ _ｉは観測時刻である（１≦ｊ≦ｑ）。前述の実施形態の各ステップにおける処理では、生の説明変数Ｘ^ｊ _ｉに代えて、各ブロックにおけるＸ^ｊ _ｉの平均値Ｚ^ｊ _ｋを用いる。このように説明変数に平均化処理を施すことにより、説明変数に混入したマイクロストラクチャー・ノイズの影響を緩和することができ、推定精度をさらに向上させることができる。

（株価の分散・共分散予測への応用）
これまで述べてきたように、日内の株価データは、高頻度非同期データの典型的な例である。本発明を用いた株価の分散・共分散の推定は、ポートフォリオのリスク・コントロール等、様々な局面に応用することができる。
例えば一般に金融機関は、保有する株式の株価変動による金融資産価値の変動をコントロールする必要がある。このため金融機関は、株価データや財務データから統計解析を用いて将来の株価の分散・共分散を予測してリスク・コントロールを行っている。特に年金基金の株式運用やインデックス・ファンドの運用においては、株価変動リスクをより正確にコントロールすることが求められるため、株価の分散・共分散を高い精度で予測することが必要となる。このような目的に対して、本発明の手法を用いた、日内株価データに基づく分散・共分散の予測は極めて有効であることが期待される。

大規模災害や金融危機などにより金融市場を取り巻く状況が急変した場合、株価のボラティリティが急激に上昇することはよく知られている。このような場合、終値等の低頻度データだけでは値動きの予測は非常に困難となる。このような場合にも、日内の株価データを観測し、本発明の手法を用いて統計解析することにより、分散・共分散の予測精度を向上させることができる。

その他にも本発明の手法は、株式市場監視の現場で、異常取引や不正取引の検知などに応用することもできる。すなわち、価格操作などが行われる傾向のあるとき、本発明の手法を用いて株価の分散・共分散を解析することにより、特定の投資家が株価ボラティリティを上げるなどの動きを検知することができる。

（電力の需要予測）
高頻度非同期データの別の例として、電力消費量のデータがある。特に猛暑期や特定のイベントがあるときなどは電力需要が大きく変わるため、これを予測することは重要である。本発明は、時々刻々モニタされる電力消費量データに基づいて、電力需要の分散・共分散を予測することに応用することができる。

（ネットワーク・サーバの需要予測）
ネットワーク・サーバの運用者にとって、サーバの需要を予測することは、サーバの効率的管理と安定的なサービス提供のために重要な課題である。本発明は、時々刻々モニタされるサーバの利用状況のデータに基づいて、サーバ需要の分散・共分散を予測することに応用することができる。

（実験例１）
図４Ａ−４Ｃおよび図５Ａ−５Ｃは、東京証券取引所における実際の日内株価データを基に、本発明の手法を用いて株価のボラティリティを予測した結果を示す。図４Ａ−４Ｃは午前の取引時間帯におけるデータ、図５Ａ−５Ｃは午後の取引時間帯におけるデータである。Ａ、ＢおよびＣの各社は、自動車関連の大型銘柄である。各図とも、株価のボラティリティの大きさを、規格化した時刻ｔおよび価格ｘの関数として予測したものである。午前の取引時間帯での予測結果からは、朝の取引開始時にボラティリティが高く、その後午前の取引終了時刻に向けて単調に減少する傾向が見られる。午後の取引時間帯での予測結果からは、やはり取引開始時にボラティリティが高く、その後単調に減少し、銘柄によっては取引終了間際に上昇する傾向が見られる。これらの傾向は、実際に市場で観測される現象と整合している。本実験例によれば、本手法を用いたボラティリティ予測が、日内株価のデータ構造を正しく分析していることが理解できる。

（実験例２）
図６は、２０１１年３月１０日の東日本震災前後における、ある電力会社株のＲＶについて、３０分間隔の実データに基づく算出値と本発明の手法を用いた予測値とを比較した結果である。実データに基づく算出値によると、３月１０日以降急激に、極めて大きなＲＶが発生していることが分かる。本手法による予測値も、実データに基づく算出値に追随して、３月１０日以降の高いＲＶを示している。従来の１日１回程度の低頻度データによる予測では、ＲＶの変動を予測するのに数ヶ月分のデータを要していた。本実験例によれば、本手法を用いたボラティリティ予測値が、ボラティリティの急激な上昇に対して迅速に反応していることが理解できる。

以上、本発明を実施例を基に説明した。これらの実施例は例示であり、それらの各構成要素や各処理プロセスの組合せにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。

高頻度非同期データには上に挙げたものの他、例えば、為替レート、放送の視聴率、インターネットサイトへのアクセス数、医療用センサからの身体データ、気象データ、地震データ、自動車の交通量、選挙の開票データ等様々なものがある。本発明による手法はこうした様々なデータに適用が可能であり、産業上の利用性が高いものである。

１０ 時系列データを取得するステップ、 ２０ 説明変数を取得するステップ、 ３０ パラメータを最適化するステップ、 ３２ 疑似対数尤度関数のバイアスを補正するステップ、 ３４ 高速近似計算を用いて逆行列を計算するステップ、 ４０ 時系列データの分散・共分散を計算するステップ

Claims (10)

1. 時系列データを取得するステップと、
   前記時系列データの時間変化の要因である説明変数を取得するステップと、
   前記時系列データと前記説明変数とに基づいて、疑似対数尤度関数を最大化することにより、前記時系列データの分散を決めるパラメータを最適化するステップと、
   最適化された前記パラメータに基づいて、前記時系列データの分散を推定するステップと、
   を備える、データ処理方法。
2. 前記パラメータを最適化するステップは、バイアスを補正するステップをさらに含む
   ことを特徴とする、請求項１に記載のデータ処理方法。
3. 前記パラメータを最適化するステップは、高速近似計算を用いて逆行列を計算するステップをさらに含む
   ことを特徴とする、請求項１または２に記載のデータ処理方法。
4. 前記時系列データは複数種類の時系列データであり、
   前記時系列データの分散を推定するステップは、前記時系列データの共分散を推定するステップをさらに含む
   ことを特徴とする、請求項１乃至３のいずれか一項に記載のデータ処理方法。
5. 前記パラメータを最適化するステップは、ニューラル・ネットワークを用いる
   ことを特徴とする、請求項１乃至４のいずれか一項に記載のデータ処理方法。
6. 前記説明変数を平均化するステップをさらに備える、
   請求項１乃至５のいずれか一項に記載のデータ処理方法。
7. 前記パラメータを最適化するステップは、高速近似計算を用いて、行列式の対数を計算するステップをさらに含む
   ことを特徴とする、請求項１乃至６のいずれか一項に記載のデータ処理方法。
8. 前記時系列データは、株価データである
   ことを特徴とする、請求項１乃至７のいずれか一項に記載のデータ処理方法。
9. 時系列データを取得する時系列データ取得部と、
   前記時系列データの時間変化の要因である説明変数を取得する説明変数取得部と、
   前記時系列データと前記説明変数とに基づいて、疑似対数尤度関数を最大化することにより、時系列データの分散を決めるパラメータを最適化するパラメータ最適化部と、
   最適化された前記パラメータに基づいて、前記時系列データの分散を推定する分散推定部と、
   を備える、データ処理装置。
10. 時系列データを取得するステップと、
    前記時系列データの時間変化の要因である説明変数を取得するステップと、
    前記時系列データと前記説明変数とに基づいて、疑似対数尤度関数を最大化することにより、前記時系列データの分散を決めるパラメータを最適化するステップと、
    最適化された前記パラメータに基づいて、前記時系列データの分散を推定するステップと、
    をコンピュータに実行させる、プログラム。