JP2014041565A - 時系列データ解析装置、方法、及びプログラム - Google Patents

Abstract

【課題】突発的な変動も考慮した時系列データの解析を行う。
【解決手段】状態空間モデリング部３１が、突発変動成分を含む複数の変動成分からなる状態ベクトルＸをシステムノイズを用いて非線形に時間更新するための状態更新式、及び状態ベクトルＸと観測値との関係を示す観測方程式によって定められた状態空間モデルであって、状態空間モデルに基づいて得られる観測値の予測値と実際の観測値との予測誤差が大きいほど大きく切り替わる突発変動ノイズを用いた時間減衰過程により突発変動成分をモデル化した状態空間モデルを定義し、状態推定部３２が、状態空間モデル及び観測値の時系列データに基づいて、状態ベクトルＸを逐次推定し、信号抽出部３４が、推定された状態ベクトルＸから、複数の変動成分の各々を示す信号を抽出する。
【選択図】図１

Description

本発明は、時系列データ解析装置、方法、及びプログラムに関する。

現実に観測される時系列データは、変動のパターンが時間と共に大きく変化する非定常性を有している。従来では、この非定常性を、主に長期変動（トレンド）成分、及び長周期の周期変動成分で捉えられてきた（例えば、非特許文献１参照）。非特許文献１では時系列データの長期変動（トレンド）を精度良く解析するため、m^th-order stochastic trendと呼ばれるモデルを提案している。このモデルでは、ｍ＝１で通常のrandom walkモデルに一致し、ｍ≧２ではトレンドの傾き成分をrandom walkでモデル化している。トレンドの推定は周波数領域で行い、最小分散推定量としてWiener-Kolmogorov filterを導出している。

Andrew C. Harvey, et al., "General Model-Based Filters for Extracting Cycles and Trends in Economic Time Series", The Review of Economics and Statistics, 2003.

しかしながら、従来手法では、突発的に生じる変動に対して、妥当な解析結果を得ることができない場合がある、という問題があった。

本発明は、上記の事情を鑑みてなされたもので、突発的な変動も考慮した時系列データの解析を行うことができる時系列データ解析装置、方法、及びプログラムを提供することを目的とする。

上記目的を達成するために、本発明の時系列データ解析装置は、突発変動成分を含む複数の変動成分からなる状態ベクトルＸをシステムノイズを用いて非線形に時間更新するための状態更新式、及び状態ベクトルＸと観測値との関係を示す観測方程式によって定められた状態空間モデルであって、前記状態空間モデルに基づいて得られる前記観測値の予測値と実際の観測値との予測誤差が大きいほど大きくなる突発変動ノイズを用いた時間減衰過程により前記突発変動成分をモデル化した状態空間モデルと、前記観測値の時系列データとに基づいて、前記状態ベクトルＸを逐次推定する状態推定手段と、前記状態推定手段により推定された状態ベクトルＸから、前記複数の変動成分の各々を示す信号を抽出する信号抽出手段と、を含んで構成されている。

本発明の時系列データ解析装置によれば、状態空間モデルとして、突発変動成分を含む複数の変動成分からなる状態ベクトルＸをシステムノイズを用いて非線形に時間更新するための状態更新式、及び状態ベクトルＸと観測値との関係を示す観測方程式によって定められたモデルを定義する。この状態ベクトルＸに含まれる突発変動成分を、状態空間モデルに基づいて得られる観測値の予測値と実際の観測値との予測誤差が大きいほど大きくなる突発変動ノイズを用いた時間減衰過程によりモデル化する。そして、状態推定手段が、この状態空間モデルと観測値の時系列データとに基づいて、状態ベクトルＸを逐次推定し、信号抽出手段が、状態推定手段により推定された状態ベクトルＸから、複数の変動成分の各々を示す信号を抽出する。

このように、時系列データの突発的な変動を示す突発変動成分を、予測誤差が大きいほど大きくなる突発変動ノイズを用いた時間減衰過程によりモデル化した状態空間モデルに基づいて、突発変動成分を含む状態ベクトルを推定するため、突発的な変動も考慮した時系列データの解析を行うことができる。

また、前記状態推定手段は、前記突発変動ノイズを、前記予測誤差をパラメータとするシグモイド関数で表すことができる。シグモイド関数は、信号の切替えによく用いられる関数であり、突発的な変動をモデル化するために適した関数である。

また、前記状態推定手段は、前記状態空間モデル及び前記観測値の時系列データに基づいて、シグマ点列の発生、状態ベクトルの状態更新、及び観測更新を行うＵＦＫ（ＵｎｓｃｅｎｔｅｄＫａｌｍａｎＦｉｌｔｅｒ）によって、前記状態ベクトルＸ及び該状態ベクトルの推定誤差分散の共分散行列Ｐを逐次推定する際に、前記状態更新後の前記共分散行列Ｐの予測値￣Ｐをシステムノイズの共分散行列Ｑの対角成分の総和で正規化した正規化共分散行列￣Ｐ’、前記観測更新後の前記共分散行列Ｐの推定値＾Ｐを前記総和で正規化した正規化共分散行列＾Ｐ’、及び前記システムノイズの共分散行列Ｑを前記総和で正規化した正規化システムノイズＱ’を用いると共に、前記総和を用いて変更された重みに基づいて、前記シグマ点列を発生させることができる。これにより、ノイズの統計量が未知であっても、通常のＵＦＫと同程度の精度で、ロバストな解析結果を得ることができる。

また、前記状態推定手段の推定結果に基づいて、前記正規化システムノイズＱ’の各分散値の比率、及び前記突発変動ノイズに関するパラメータをパラメトライズして得られる尤度が最大となるときのパラメータを、前記状態空間モデルのパラメータとして同定するモデル同定手段を含んで構成することができる。

また、前記状態推定手段は、前記総和が収束するまで前記ＵＦＫを繰り返し行い、前記総和が収束した際の前記状態ベクトルＸの推定値＾Ｘを平滑化した平滑化推定値~Ｘを推定することができる。

また、本発明の時系列データ解析方法は、状態推定手段と、信号抽出手段とを含む時系列データ解析装置における時系列データ解析方法であって、前記状態推定手段が、突発変動成分を含む複数の変動成分からなる状態ベクトルＸをシステムノイズを用いて非線形に時間更新するための状態更新式、及び状態ベクトルＸと観測値との関係を示す観測方程式によって定められた状態空間モデルであって、前記状態空間モデルに基づいて得られる前記観測値の予測値と実際の観測値との予測誤差が大きいほど大きくなる突発変動ノイズを用いた時間減衰過程により前記突発変動成分をモデル化した状態空間モデルと、前記観測値の時系列データとに基づいて、前記状態ベクトルＸを逐次推定し、前記信号抽出手段が、前記状態推定手段により推定された状態ベクトルＸから、前記複数の変動成分の各々を示す信号を抽出する方法である。

また、本発明の時系列データ解析プログラムは、コンピュータを、上記の時系列データ解析装置を構成する各手段として機能させるためのプログラムである。

以上説明したように、本発明の時系列データ解析装置、方法、及びプログラムによれば、時系列データの突発的な変動を示す突発変動成分を、予測誤差が大きいほど大きくなる突発変動ノイズを用いた時間減衰過程によりモデル化した状態空間モデルに基づいて、突発変動成分を含む状態ベクトルを推定するため、突発的な変動も考慮した時系列データの解析を行うことができる、という効果が得られる。

本実施の形態に係る時系列データ解析装置の構成を示す概略図である。 観測値の時系列データのデータ構造の一例を示す図である。 出力データのデータ構造の一例を示す図である。 本実施の形態における時系列データ解析処理ルーチンの内容を示すフローチャートである。 実験に用いた時系列データを示すグラフである。 本実施の形態のモデルによる循環変動成分の抽出結果を示すグラフである。 本実施の形態のモデルによる長期変動成分の抽出結果を示すグラフである。 本実施の形態のモデルによる周期変動成分の抽出結果を示すグラフである。 本実施の形態のモデルによる突発変動成分の抽出結果を示すグラフである。 既存モデルによる循環変動成分の抽出結果を示すグラフである。 本実施の形態のモデルの突発変動ノイズとして観測ノイズを適用した場合のモデルによる循環変動成分の抽出結果を示すグラフである。

以下、図面を参照して本発明の実施の形態を詳細に説明する。

＜本実施の形態の概要＞
本実施の形態に係る時系列データ解析装置では、まず、時系列データの突発的な変動を時間減衰過程としてモデル化する。そして通常の時系列成分モデルにこの突発変動成分を加えたモデルを状態空間で表現する。状態の推定は、後述するＮｏｉｓｅｓ Ｎｏｒｍａｌｉｚｅｄ Ｕｎｓｃｅｎｔｅｄ Ｋａｌｍａｎ Ｆｉｌｔｅｒ ＆ ＵＲＴＳ Ｓｍｏｏｔｈｅｒ（ＮＮ−ＵＫＦ＆ＵＲＴＳＳ）を使用して行い、突発的な事象をオンラインで検出しながら推定を実行する。

＜システム構成＞
本実施の形態に係る時系列データ解析装置１０は、ＣＰＵと、ＲＡＭと、後述する時系列データ解析処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成されている。このコンピュータは、機能的には、図１に示すように、入力部２０と、演算部３０と、出力部４０とを含んだ構成で表すことができる。

入力部２０は、入力された観測値の時系列データを受け付ける。観測値は、スカラー値であり、すなわち、時系列データは、１次元時系列である。観測値は、例えば、図２に示すように、時刻ｔと観測値ｙ_ｔとの組（ｔ＝１，・・・，Ｎ）からなる時系列データである。

また、入力部２０は、オペレータから、後述する状態空間モデルを定義するための設定値の入力を受け付ける。

演算部３０は、状態空間モデリング部３１と、状態推定部３２と、モデル同定部３３と、信号抽出部３４とを含んだ構成で表すことができる。

状態空間モデリング部３１は、オペレータにより入力部２０から入力された設定値を受け付けて、下記（１）式及び（２）式に示す状態空間モデルを定義する。

ただし、Ｘ_ｔは時刻ｔにおける状態ベクトルで、通常、非観測量が取られる。状態ベクトルＸ_ｔを時系列データｙ_ｔを用いて推定することを目的とする。μはシステムノイズベクトルであり、正規分布に従う白色ノイズベクトルである。また、Ｎ（ｍ、ｖ）は平均ｍ、分散ｖの正規分布を表す。したがって、Ｑはシステムノイズの共分散行列に相当する。なお、ここでは、観測ノイズフリーとしている。上記（２）式に示されるように、ｙ_ｔは時刻ｔにおける観測値で、状態Ｘ_ｔの線形変換Ｈにより生成される。上記（２）式は、状態ベクトルＸ_ｔと観測値ｙ_ｔとの関係を示している。また、上記（１）式のｆ（・）が非線形箇所であり、これによりシステムモデルは非線形になる。ｆ（・）の型はモデルにより異なる。

次に、突発的な変動Ｅは、下記（３）式に示す時間減衰過程としてモデル化する。

ただし、τ_ｔは時定数であり、Ｅ_ｔが１／ｅになるまでの時間を意味する。また、μ'はシステムノイズμの共分散行列の対角成分、すなわちシステムノイズμの各分散値にＧ_ｔを掛けたものである。（３）式の離散解は、下記（４）式となる。

ただし、Δは時系列データの観測幅であり、本実施の形態では１として説明する。この突発変動成分を加えた時系列成分モデルの観測方程式は、下記（５）式となる。

ただし、ｘ_ｔは時刻ｔにおける循環変動成分（ＡＲ過程）、Ｔ_ｔは長期変動（トレンド）成分、β_ｔは周期変動成分、Ｅ_ｔは突発変動成分である。それぞれの成分の時間更新式は、下記（６）式〜（１０）式となる。

循環変動成分はＭ次のＡＲ過程、長期変動（トレンド）成分は線形トレンド、周期変動成分は周期Ｃでモデル化してある。時定数τ_ｔは突発変動成分Ｅ_ｔをある程度の自己相関を持った減衰過程として抽出するため、例えば、τ_ｔ∈［７，１３］という制限を持った局所一定モデルに従い時間発展するものとする。なお、μ_ｉ（ｉ＝１，２，・・・，２Ｍ＋３＋Ｃ＋１）は、システムノイズμのｉ番目の分散値であることを表す。
状態ベクトルの定義は下記（１１）式となり、各変動成分と関連パラメータ（ＡＲ係数及び時定数）とを同時に推定する。

観測モデル（２）式のＨはＨ＝［1,0,・・・,1,0,1,0,・・・,0,1,0］であり、（５）式の観測方程式と同等になるように設計する。

また、システムモデル（１）式のＧ_ｔは、下記（１２）式となる。

ここで、ν_ｔは時刻ｔにおける予測誤差で、この予測誤差が非常に大きくなった時点を突発的な変動の開始とみなし、ｇ（ν_ｔ）が大きくなる様に設定する。その結果、（９）式のμ'_{（２Ｍ＋３＋Ｃ）}＝ｇ（ν_ｔ）μ_{（２Ｍ＋３＋Ｃ）}が大きくなり、Ｅ_ｔが突発変動を成分として推定し始めることになる。

状態推定部３２は、入力された時系列データと、状態空間モデリング部３１でモデリングされた状態空間モデルとに基づき、状態ベクトルの逐次推定を行う。モデルは非線形であるので、以下に示すＮＮ−ＵＫＦ＆ＵＲＴＳＳによって、状態ベクトルＸ_ｔの逐次状態推定（オンライン状態推定）を行う。ＮＮ−ＵＫＦのアルゴリズムは、下記（１３）式〜（２５）式で表されるシグマ点列の発生、状態更新（１期先予測）、及び観測更新（フィルタ）からなる。

ただし、κは予め定められた定数であり、ｎは状態ベクトルＸ_ｔの次元数である。Ｐ_ｔは状態ベクトルＸ_ｔの推定誤差分散の共分散行列である。また、￣Ａ（数式中では、記号Ａ上に“￣”）は状態更新により予測されたＡの予測値を示し、＾Ａ（数式中では、記号Ａ上に“＾”）は観測更新により推定されたＡの推定値を示す。なお、￣Ｐ'_ｔ、＾Ｐ'_ｔ、Ｑ’は下記（２６）式〜（２８）式に示すように、共分散行列の予測値￣Ｐ_ｔ、共分散行列の推定値＾Ｐ_ｔ、及びシステムノイズの共分散行列Ｑを、ｔｒ（Ｑ）で割って正規化したものである。なお、観測方程式に観測ノイズを含む場合には、ＮＮ−ＵＫＦアルゴリズムに観測ノイズの分散値に関するｒ’の項が含まれるが、本実施の形態では、観測ノイズフリーとしているため、ＮＮ−ＵＫＦアルゴリズム内のｒ’は０としている。

ここで、ｔｒ（Ｑ）はシステムノイズμの共分散行列Ｑの対角成分、すなわちシステムノイズμの分散値の総和であり、アルゴリズム中のαとの関係はｔｒ（Ｑ）＝１／α^２である。多くの場合、ｔｒ（Ｑ）の値は未知であるので、最初はｔｒ（Ｑ）＝１、すなわちα＝１に設定してＮＮ−ＵＫＦを実行する。そして、ＮＮ−ＵＫＦの実行毎に、下記（３０）式及び（３１）式に基づいて、ｔｒ（Ｑ）を推定する。

状態推定部３２は、ＮＮ−ＵＫＦ実行時に使用したｔｒ（Ｑ）_ｎｏｗの値と、ＮＮ−ＵＫＦ実行終了時に推定したｔｒ（Ｑ）_ｎｅｗの値との差が所定の閾値以下であれば、ｔｒ（Ｑ）の値が収束したとみなし、この実行によって推定された状態の結果を採用する手法を取る。

また、（１２）式のｇ（ν_ｔ）は、下記（３２）式〜（３４）式となる。

ここで、σ_ｔは時刻ｔにおける予測分布の標準偏差で、ｚ_ｔは予測誤差ν_ｔの大きさを表す量である。ξはｚ_ｔの平均を表すパラメータであり、本実施の形態では、ξ＝４，５と振っている。ｇ（ν_ｔ）は標準シグモイド関数で、０〜１の値を取る。シグモイド関数はその形から信号の切替えによく使われ、本実施の形態においても、上述したように、突発変動成分の時間発展時にのるノイズ（突発変動ノイズ）の量の切替えに使用する。通常の予測誤差の範囲であればｇ（ν_ｔ）は０に近い値を取り、従って、μ'_{（２Ｍ＋３＋Ｃ）}＝ｇ（ν_ｔ）μ_{（２Ｍ＋３＋Ｃ）}も小さくなり、突発変動成分Ｅ_ｔは正常に減衰する。一方、予測誤差が大きく外れる様なデータを観測すると、ｇ（ν_ｔ）は１に近い値を取り、μ'_{（２Ｍ＋３＋Ｃ）}＝ｇ（ν_ｔ）μ_{（２Ｍ＋３＋Ｃ）}も大きくなり、突発変動成分Ｅ_ｔは観測更新（フィルタ）によりその値が大きく修正されることになる。この予測誤差に基づくノイズの切替えにより、突発変動成分を自己相関減衰過程として抽出する。

最後の観測値ｙ_ｔ＝Ｎまでの処理が完了すると、（２１）式を使用したＮｏｉｓｅｓ Ｎｏｒｍａｌｉｚｅｄ Ｕｎｓｃｅｎｔｅｄ ＲＴＳ Ｓｍｏｏｔｈｅｒ（ＮＮ−ＵＲＴＳＳ）により、下記（３５）式〜（３７）式に示す平滑化を行う。

ここで、Ｓ_ｔ−１は平滑化ゲイン、~Ｘ_ｔ−１及び~Ｐ_ｔ−１（数式中では、記号ＸまたはＰ上に“~”）はそれぞれ状態の期待値及び共分散行列の平滑化推定値である。

モデル同定部３３は、状態推定部３２により正規化されたシステムノイズの共分散行列Ｑ’の対角成分（分散値）の総和μ_１＋μ_２＋・・・＋μ_{（２Ｍ＋３＋Ｃ＋１）}を利用して、各分散値の比率γ、及び（３３）で導入したξ（例えば、ξ＝４，５）をパラメトライズして振って、γとξとの組（γ，ξ）毎に下記（３８）式〜（４１）式に示す対数尤度ＬＬ（γ，ξ）を求める。なお、（６）式及び（９）式に示した循環変動成分及び突発変動成分にのるノイズμ_１及びμ'_{（２Ｍ＋３＋Ｃ）}以外は非常に小さい（?０）ため、ここでは、γについては、μ_１とμ'_{（２Ｍ＋３＋Ｃ）}との比率のみをパラメトライズする。そして、対数尤度ＬＬ（γ、ξ）が最大になった（γ，ξ）から得られるＱ’を、状態推定モデルのパラメータとして推定することにより、モデルの同定を行う。

ここで、Ｍ（γ，ξ）がシステムノイズμの分散値の総和の推定値となる。

信号抽出部３４は、状態推定部３２により推定された平滑化推定値~Ｘ_ｔ−１から、各成分（循環変動成分ｘ_ｔ、長期変動（トレンド）成分Ｔ_ｔ、周期変動成分β_ｔ、突発変動成分Ｅ_ｔ）を示す信号（信号Ａ_ｔ，信号Ｂ_ｔ，・・・）を抽出する。

出力部４０は、モデル同定部３３により同定されたモデルを示すパラメータの推定値、及び信号抽出部３４により抽出された信号を出力する。信号抽出部３４により抽出された信号は、例えば、図３に示すように、時刻ｔと抽出された信号Ａ_ｔ，信号Ｂ_ｔ・・・との組（ｔ＝１，・・・，Ｎ）からなる時系列データとして出力することができる。

＜時系列データ解析装置の作用＞
次に、本実施の形態に係る時系列データ解析装置１０の作用について説明する。まず、オペレータにより、状態空間モデルを定義するための設定値が時系列データ解析装置１０に入力されると、時系列データ解析装置１０によって、入力された設定値が、メモリ（図示省略）へ格納される。また、時刻ｔ_１〜ｔ_Ｎの観測値からなる時系列データが、時系列データ解析装置１０に入力されると、時系列データ解析装置１０によって、入力された時系列データが、メモリへ格納される。そして、時系列データ解析装置１０によって、図４に示す時系列データ解析処理ルーチンが実行される。

まず、ステップ１００で、メモリに格納された状態空間モデルを定義するための設定値及び観測値の時系列データを取得する。そして、ステップ１０２で、状態空間モデリング部３１が、取得した設定値を用いて、上記（１）式及び（２）式に示す状態空間モデルを定義する。

次に、ステップ１０４で、状態推定部３２が、ＮＮ−ＵＫＦ＆ＵＲＴＳＳを実行するための初期設定として、システムノイズの分散値の総和（ｔｒ（Ｑ））、及び時系列データの時刻を示すｔをそれぞれ１にセットする。

次に、ステップ１０６で、状態推定部３２が、ＮＮ−ＵＫＦのアルゴリズムに従って、上記ステップ１００で取得した観測値の時系列データ、及び上記ステップ１０４で設定された初期設定を用いて、シグマ点列の発生、状態更新、及び観測更新を逐次行って、状態ベクトル及びその共分散行列を推定する。

次に、ステップ１０８で、状態推定部３２が、上記（３０）式及び（３１）式によりｔｒ（Ｑ）を推定し、上記ステップ１０６で使用したｔｒ（Ｑ）_ｎｏｗの値と、本ステップで推定したｔｒ（Ｑ）_ｎｅｗの値との差が所定の閾値以下か否かを判定することにより、ｔｒ（Ｑ）の値が収束したか否かを判定する。ｔｒ（Ｑ）の値が収束していない場合には、ステップ１０６へ戻り、推定された状態ベクトル＾Ｘ_ｔ、推定され正規化された共分散行列＾Ｐ'_ｔ、及び推定されたｔｒ（Ｑ）を用いて、ＮＮ−ＵＫＦによる状態ベクトル及びその共分散行列の推定を繰り返す。

一方、ｔｒ（Ｑ）の値が収束した場合には、ステップ１１０へ移行し、ｔがＮになったか否かを判定することにより、最後の観測値ｙ_Ｎまで処理が完了したか否かを判定する。ｔがＮに達していない場合には、ステップ１１２へ移行して、ｔを１インクリメントして、ステップ１０６へ戻る。ｔ＝Ｎとなった場合には、ステップ１１４へ移行する。

ステップ１１４では、状態推定部３２が、上記（２１）式を使用したＮＮ−ＵＲＴＳＳにより、上記（３５）式〜（３７）式に示す平滑化を行い、状態ベクトルの平滑化推定値~Ｘ_ｔ−１及びその共分散行列の平滑化推定値~Ｐ_ｔ−１を得る。

次に、ステップ１１６で、モデル同定部３３が、正規化されたシステムノイズの共分散行列Ｑ’の対角成分（分散値）の総和μ_１＋μ_２＋・・・＋μ_{（２Ｍ＋３＋Ｃ＋１）}を利用して、各分散値の比率γ、及び（３３）で導入したξ（ξ＝４，５）をパラメトライズして振って、γとξとの組（γ，ξ）毎に（３８）式〜（４１）式に示す対数尤度ＬＬ（γ，ξ）を求める。そして、対数尤度ＬＬ（γ、ξ）が最大になった（γ，ξ）から得られるＱ’を、状態推定モデルのパラメータとして推定することにより、モデルの同定を行う。

次に、ステップ１１８で、信号抽出部３４が、上記ステップ１１４で推定された平滑化推定値~Ｘ_ｔ−１から、各成分（循環変動成分ｘ_ｔ、長期変動（トレンド）成分Ｔ_ｔ、周期変動成分β_ｔ、突発変動成分Ｅ_ｔ）を示す信号（信号Ａ_ｔ，信号Ｂ_ｔ，・・・）を抽出する。

次に、ステップ１２０で、出力部４０が、上記ステップ１１６で同定されたモデルを示すパラメータの推定値、及び上記ステップ１１８で抽出された信号を出力し、時系列データ解析処理ルーチンを終了する。

＜実験結果＞
ここで、あるＥＣサイトの日毎の売上げ高の時系列データ（図５）を、本実施の形態における状態推定モデルで解析した実験結果を示す。図６は循環変動成分の抽出結果、図７は長期変動成分の抽出結果、図８は周期変動成分の抽出結果、図９は突発変動成分の抽出結果である。ここでは、循環変動成分のＡＲ次数Ｍは３、周期変動成分の周期Ｃは７、つまり１週間周期とした。

次に、既存のモデル（循環変動成分、長期変動成分、周期変動成分からなる状態ベクトルを推定するモデル、突発変動成分を含まないモデル）に基づき、上述の時系列データを解析した際の、循環変動成分の抽出結果を図１０に、観測ノイズの分散値を本実施の形態と同様な手法で切替えたモデル、すなわち、突発的な変動は大きな観測ノイズがのったことに依るものと見なすモデルに基づく抽出結果を図１１に示す。本実施の形態に観測ノイズを適用する場合には、（２）式の観測モデルに観測ノイズｗ〜Ｎ（０，ｒ）を加え、（３）式及び（４）式の突発変動成分の定義のμ’をｗ’（ｗ’＝Ｇ_ｔｗ）とする。また、ＮＮ−ＵＫＦアルゴリズムにおいてｒ’＝０とした箇所に適切なｒ’を用いる。

抽出された循環変動成分はそれぞれ明らかに非定常になっており、推定されたＡＲ係数も非定常と定常の境界に近い値を示していた。本実施の形態のモデルの対数尤度は-18063で、既存モデルの対数尤度はそれぞれ-19069、-18311であることから、突発的な変動を含む非定常時系列データに対して本実施の形態のモデルが有効であることを示している。また、突発的な変動を含まないデータに対しては突発変動成分Ｅ_ｔがほぼ０として推定されるため、本実施の形態のモデルは様々な時系列データに対して汎用的に使用可能である。

以上説明したように、時系列データの突発的な変動を示す突発変動成分を、予測誤差が大きいほど大きく切り替わるノイズをのせた時間減衰過程としてモデル化した状態空間モデルに基づいて、突発変動成分を含む状態ベクトルを推定するため、突発的な変動も考慮した時系列データの解析を行うことができる。

現実の時系列データは、時間と共にその様相が大きく変化する非定常性を有している。カルマンフィルタに代表されるフィルタリング手法は時系列データを生成する状態を確率的に推定し、特に時変する状態の追従性に優れている。本実施の形態では、既存の時系列モデルでは取り扱ってこなかった突発的な変動をモデル化することにより、突発事象の形相を原系列から分離抽出することが可能になり、突発変動成分のさらなる詳細な分析が可能となる。また、本実施の形態の突発変動成分を既存の成分モデルに加えることで、変化が激しい時系列データに対しても安定的な状態推定を行うことができる。

なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

例えば、状態空間モデルの定義を予め設定しておき、状態空間モデルを定義するための設定値の入力を不要としてもよい。

また、本実施の形態では、状態空間モデルのパラメータ推定値及び抽出した信号を出力する場合について説明したが、状態推定部により推定された平滑化推定値~Ｘ_ｔ−１及び~Ｐ_ｔ−１を出力するようにしてもよいし、￣ｙ＝Ｈ￣Ｘに従って観測値の予測値を求めて出力するようにしてもよい。

また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。

１０ 時系列データ解析装置
２０ 入力部
３０ 演算部
３１ 状態空間モデリング部
３２ 状態推定部
３３ モデル同定部
３４ 信号抽出部
４０ 出力部

Claims (7)

1. 突発変動成分を含む複数の変動成分からなる状態ベクトルＸをシステムノイズを用いて非線形に時間更新するための状態更新式、及び状態ベクトルＸと観測値との関係を示す観測方程式によって定められた状態空間モデルであって、前記状態空間モデルに基づいて得られる前記観測値の予測値と実際の観測値との予測誤差が大きいほど大きくなる突発変動ノイズを用いた時間減衰過程により前記突発変動成分をモデル化した状態空間モデルと、前記観測値の時系列データとに基づいて、前記状態ベクトルＸを逐次推定する状態推定手段と、
   前記状態推定手段により推定された状態ベクトルＸから、前記複数の変動成分の各々を示す信号を抽出する信号抽出手段と、
   を含む時系列データ解析装置。
2. 前記状態推定手段は、前記突発変動ノイズを、前記予測誤差をパラメータとするシグモイド関数で表した請求項１記載の時系列データ解析装置。
3. 前記状態推定手段は、前記状態空間モデル及び前記観測値の時系列データに基づいて、シグマ点列の発生、状態ベクトルの状態更新、及び観測更新を行うＵＦＫ（Ｕｎｓｃｅｎｔｅｄ Ｋａｌｍａｎ Ｆｉｌｔｅｒ）によって、前記状態ベクトルＸ及び該状態ベクトルの推定誤差分散の共分散行列Ｐを逐次推定する際に、前記状態更新後の前記共分散行列Ｐの予測値￣Ｐをシステムノイズの共分散行列Ｑの対角成分の総和で正規化した正規化共分散行列￣Ｐ’、前記観測更新後の前記共分散行列Ｐの推定値＾Ｐを前記総和で正規化した正規化共分散行列＾Ｐ’、及び前記システムノイズの共分散行列Ｑを前記総和で正規化した正規化システムノイズＱ’を用いると共に、前記総和を用いて変更された重みに基づいて、前記シグマ点列を発生させる請求項１または請求項２記載の時系列データ解析装置。
4. 前記状態推定手段の推定結果に基づいて、前記正規化システムノイズＱ’の各分散値の比率、及び前記突発変動ノイズに関するパラメータをパラメトライズして得られる尤度が最大となるときのパラメータを、前記状態空間モデルのパラメータとして同定するモデル同定手段を含む請求項１〜請求項３のいずれか１項記載の時系列データ解析装置。
5. 前記状態推定手段は、前記総和が収束するまで前記ＵＦＫを繰り返し行い、前記総和が収束した際の前記状態ベクトルＸの推定値＾Ｘを平滑化した平滑化推定値~Ｘを推定する請求項１〜請求項４のいずれか１項記載の時系列データ解析装置。
6. 状態推定手段と、信号抽出手段とを含む時系列データ解析装置における時系列データ解析方法であって、
   前記状態推定手段が、突発変動成分を含む複数の変動成分からなる状態ベクトルＸをシステムノイズを用いて非線形に時間更新するための状態更新式、及び状態ベクトルＸと観測値との関係を示す観測方程式によって定められた状態空間モデルであって、前記状態空間モデルに基づいて得られる前記観測値の予測値と実際の観測値との予測誤差が大きいほど大きくなる突発変動ノイズを用いた時間減衰過程により前記突発変動成分をモデル化した状態空間モデルと、前記観測値の時系列データとに基づいて、前記状態ベクトルＸを逐次推定し、
   前記信号抽出手段が、前記状態推定手段により推定された状態ベクトルＸから、前記複数の変動成分の各々を示す信号を抽出する
   時系列データ解析方法。
7. コンピュータを、請求項１〜請求項５のいずれか１項記載の時系列データ解析装置を構成する各手段として機能させるための時系列データ解析プログラム。

Images