(57)【要約】【課題】複数のデータを含む時系列からより正確にパワースペクトル密度を推定する方法、パワースペクトル密度推定装置、コンピュータプログラム、及び記録媒体を提供する。【解決手段】複数のデータを含む時系列｛ｙｔ　｝から、自己相関関数の推定値｛Ｒｋ　：ｋ＝０，１，…，Ｋ｝をＲｋ　＝Ｅ［ｙｔ　ｙｔ＋ｋ　］と計算し（Ｓ２）、前記時系列にＡＲモデルを適用したときのパラメータの推定値｛αｍ　：ｍ＝１，２，…，Ｍ｝を、拡張ユール・ウォーカー方程式を用いて導かれる式により計算し（Ｓ３〜Ｓ１０）、ユール・ウォーカー方程式及び拡張ユール・ウォーカー方程式を用いて導かれる式により、前記時系列のパワースペクトル密度の近似式Ｓ（ω）を計算する（Ｓ１１〜Ｓ１２）。【選択図】　　　　図２
2004118355.tif 000002
【特許請求の範囲】
【請求項１】
記憶部及び演算部を備えたコンピュータを用いて、複数のデータを含む時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定する方法において、
時系列をなす複数のデータを記憶部に記憶し、
拡張ユール・ウォーカー方程式を用いて導かれる式により、記憶している時系列からＡＲモデルのパラメータの推定値を演算部にて計算し、
計算した前記パラメータの推定値を用いて、ユール・ウォーカー方程式及び拡張ユール・ウォーカー方程式を用いて導かれる式により、前記時系列のパワースペクトル密度の推定値を演算部にて計算する
ことを特徴とするパワースペクトル密度推定方法。
【請求項２】
解析対象から得られる時系列のデータを用いて前記解析対象の状態を解析するために、記憶部及び演算部を備えたコンピュータを用いて、前記時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定する方法において、
解析対象から得られる時系列のデータ、及び該データが得られるときのサンプリングの時間間隔を記憶部に記憶し、
拡張ユール・ウォーカー方程式を用いて導かれる式に、記憶している時系列のデータを代入して、ＡＲモデルのパラメータの推定値を演算部にて計算し、
ユール・ウォーカー方程式及び拡張ユール・ウォーカー方程式を用いて導かれる式に、計算した前記パラメータの推定値、前記時系列のデータ、及び記憶している前記時間間隔を代入して、前記時系列のパワースペクトル密度の推定値を演算部にて計算する
ことを特徴とするパワースペクトル密度推定方法。
【請求項３】
記憶部及び演算部を備えたコンピュータを用いて、複数のデータを含む時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定する方法において、
時系列をなす複数のデータ｛ｙｔ　｝、及び該データ間の時間間隔Δを記憶部に記憶し、
記憶している時系列の自己相関関数の推定値｛Ｒｋ　：ｋ＝０，１，…，Ｋ｝を、
Ｒｋ　＝Ｅ［ｙｔ　ｙｔ＋ｋ　］
と演算部にて計算し、
前記時系列にＡＲモデルを適用したときのパラメータの推定値｛αｍ　：ｍ＝１，２，…，Ｍ｝を、
【数１】
2004118355.tif 000028
と演算部にて計算し、
記憶している前記時間間隔Δ及び計算したパラメータの推定値｛αｍ　｝を用いて、前記時系列のパワースペクトル密度の角周波数ωについての近似式Ｓ（ω）を、
【数２】
2004118355.tif 000029
の式を用いて演算部にて計算する
ことを特徴とするパワースペクトル密度推定方法。
【請求項４】
記憶部及び演算部を備えたコンピュータを用いて、複数のデータを含む時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定する方法において、
時系列をなす複数のデータ｛ｙｔ　｝、及び該データ間の時間間隔Δを記憶部に記憶し、
記憶している時系列の自己相関関数の推定値｛Ｒｋ　：ｋ＝０，１，…，Ｋ｝を、
Ｒｋ　＝Ｅ［ｙｔ　ｙｔ＋ｋ　］
と演算部にて計算し、
前記時系列にＡＲモデルを適用したときのパラメータの推定値｛αｍ　：ｍ＝１，２，…，Ｍ｝を、
【数３】
2004118355.tif 000030
の式を用いて演算部にて計算し、
記憶している前記時間間隔Δ及び計算したパラメータの推定値｛αｍ　｝を用いて、前記時系列のパワースペクトル密度の角周波数ωについての近似式Ｓ（ω）を、
【数４】
2004118355.tif 000031
と演算部にて計算する
ことを特徴とするパワースペクトル密度推定方法。
【請求項５】
複数のデータを含む時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定する装置において、
時系列をなす複数のデータ｛ｙｔ　｝、及び該データ間の時間間隔Δを記憶する手段と、
記憶している時系列の自己相関関数の推定値｛Ｒｋ　：ｋ＝０，１，…，Ｋ｝を、
Ｒｋ　＝Ｅ［ｙｔ　ｙｔ＋ｋ　］
と計算する手段と、
前記時系列にＡＲモデルを適用したときのパラメータの推定値｛αｍ　：ｍ＝１，２，…，Ｍ｝を、
【数５】
2004118355.tif 000032
の式を用いて計算する手段と、
記憶している前記時間間隔Δ及び計算したパラメータの推定値｛αｍ　｝を用いて、前記時系列のパワースペクトル密度の角周波数ωについての近似式Ｓ（ω）を、
【数６】
2004118355.tif 000033
と計算する手段と
を備えることを特徴とするパワースペクトル密度推定装置。
【請求項６】
コンピュータに、複数のデータを含む時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定させるコンピュータプログラムにおいて、
コンピュータに、記憶している時系列をなす複数のデータ｛ｙｔ　｝の自己相関関数の推定値｛Ｒｋ　：ｋ＝０，１，…，Ｋ｝を、
Ｒｋ　＝Ｅ［ｙｔ　ｙｔ＋ｋ　］
と計算させる手順と、
コンピュータに、前記時系列にＡＲモデルを適用したときのパラメータの推定値｛αｍ　：ｍ＝１，２，…，Ｍ｝を、
【数７】
2004118355.tif 000034
の式を用いて計算させる手順と、
コンピュータに、記憶している前記データ間の時間間隔Δ、及び計算したパラメータの推定値｛αｍ　｝を用いて、前記時系列のパワースペクトル密度の角周波数ωについての近似式Ｓ（ω）を、
【数８】
2004118355.tif 000035
と計算させる手順と
を含むことを特徴とするコンピュータプログラム。
【請求項７】
コンピュータに、複数のデータを含む時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定させるコンピュータプログラムを記録してあるコンピュータでの読み取りが可能な記録媒体において、
コンピュータに、記憶している時系列をなす複数のデータ｛ｙｔ　｝の自己相関関数の推定値｛Ｒｋ　：ｋ＝０，１，…，Ｋ｝を、
Ｒｋ　＝Ｅ［ｙｔ　ｙｔ＋ｋ　］
と計算させる手順と、
コンピュータに、前記時系列にＡＲモデルを適用したときのパラメータの推定値｛αｍ　：ｍ＝１，２，…，Ｍ｝を、
【数９】
2004118355.tif 000036
の式を用いて計算させる手順と、
コンピュータに、記憶している前記データ間の時間間隔Δ、及び計算したパラメータの推定値｛αｍ　｝を用いて、前記時系列のパワースペクトル密度の角周波数ωについての近似式Ｓ（ω）を、
【数１０】
2004118355.tif 000037
と計算させる手順と
を含むコンピュータプログラムを記録してあることを特徴とするコンピュータでの読み取りが可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、複数のデータを含む時系列から、該時系列のパワースペクトル密度を推定する方法、パワースペクトル密度推定装置、コンピュータを該パワースペクトル密度推定装置として実現するためのコンピュータプログラム、及び記録媒体に関する。
【０００２】
【従来の技術】
離散的な複数のデータを含む時系列について、該時系列の解析、該時系列の加工、又は新たなデータの予測などを行うために、時系列の各波長または周波数毎の挙動を調べるフーリエ解析がしばしば行われる。フーリエ解析のために時系列からパワースペクトル密度を得る方法には、標準法と高速フーリエ変換（ＦＦＴ）法とがある。
【０００３】
定常エルゴード的不規則過程の一つの標本関数ｙ（ｔ）が与えられたときには、この過程の自己相関関数Ｒ（τ）及びパワースペクトル密度Ｓ（ω）は、
【０００４】
【数１１】
2004118355.tif 000003
【０００５】
で与えられる。ここでｉは虚数であり、ωは角周波数である。τ＝０，Δ，２Δ，…におけるＲ（τ）のみが得られているときには、（２）式の近似式として、
【０００６】
【数１２】
2004118355.tif 000004
【０００７】
が計算される。ｎの数をＭでうち切ることによって、（３）式は、更に、
【０００８】
【数１３】
2004118355.tif 000005
【０００９】
と近似される。ここで、Ｗｎ　はラグウインドウであり、ラグウインドウＷ（ｔ）について、Ｗ（ｎΔ）＝Ｗｎ　とした。更に、ｙｎ　＝ｙ（ｎΔ）である標本時系列ｙｎ　について、Ｒ（ｎΔ）の推定値Ｒｋ　を、
【００１０】
【数１４】
2004118355.tif 000006
【００１１】
とする。標準法では、（５）式で計算されるＲｋ　を用いて、パワースペクトル密度Ｓ（ω）を
【００１２】
【数１５】
2004118355.tif 000007
【００１３】
と近似する。（６）式を用いることでコンピュータでのパワースペクトル密度の計算が可能となる。しかし、（４）式の近似の段階で計算を有限項で打ち切ることによる偏り誤差が生じ、（６）式の近似では統計的誤差が導入される。ＦＦＴ法は、計算原理は標準法とは全く異なり、計算時間の面で標準法よりもやや優れているが、偏り誤差と統計的誤差とが生じるのは同様である。
【００１４】
パワースペクトル密度を推定する他の方法として、パワースペクトル密度を推定しようとする時系列を、ＡＲ（Ａｕｔｏ−Ｒｅｇｒｅｓｓｉｖｅ　）モデルに従った時系列であるとして推定を行う方法がある。複数のデータを含む時系列｛ｙｔ　：ｔ＝…，−１，０，１，２，…｝を、定常正規エルゴード過程の標本時系列としたときには、ｙｔ　は、
ｙｔ　＋α１　ｙｔ−１　＋…＋αＭ　ｙｔ−Ｍ　＝ｕｔ　　　…（７）
と表される。ここで、｛αｍ　：ｍ＝１，２，…，Ｍ｝はＡＲモデルのパラメータであり、｛ｕｔ　｝は平均値０及び分散σ２　の正規性雑音とする。（７）式で示される如きＡＲモデルでのスペクトル密度Ｓ（ω）の推定値は、
【００１５】
【数１６】
2004118355.tif 000008
【００１６】
で与えられる。（８）式は（３）式と等価であるため、（８）式を用いたパワースペクトル密度の推定は、（６）式を用いる標準法、又はＦＦＴ法に比べて、推定精度の点で優れている。しかし、実際に得られる情報は時系列｛ｙｔ　｝のデータのみであり、（８）式を用いたパワースペクトル密度の推定を行うためには、ＡＲモデルのパラメータ｛αｍ　｝を定める必要があり、前記パラメータの推定の際に誤差が導入される。非特許文献１では、時系列｛ｙｔ　｝に含まれるデータの値ｙｔ　とＡＲモデルによって計算される前記データの値ｙｔ　の予測値との２乗和が最小となるようにＡＲモデルのパラメータを定め、定めたパラメータを用いて前記時系列のパワースペクトル密度を推定するＦＰＥ（Ｆｉｎａｌ　Ｐｒｅｄｉｃｔｉｏｎ　Ｅｒｒｏｒ）法が開示されている。また、非特許文献２では、ＡＲＭＡ（Ａｕｔｏｒｅｇｒｅｓｓｉｖｅｍｏｖｉｎｇ　ａｖｅｒａｇｅ　）モデルに従う時系列について、これと等価なＡＲモデルを考えることによって前記時系列のパワースペクトルを推定するＦＯＥ（ｆｕｎｃｔｉｏｎ　ｏｆ　ｏｒｄｅｒ　ｅｓｔｉｍａｔｉｏｎ）法が開示されている。
【００１７】
【非特許文献１】
赤池（Ｈ．　Ａｋａｉｋｅ），「フィッティング　オートレグレッシブ　モデルズ　フォア　プレディクション（Ｆｉｔｔｉｎｇ　Ａｕｔｏｒｅｇｒｅｓｓｉｖｅ　Ｍｏｄｅｌｓ　ｆｏｒ　Ｐｒｅｄｉｃｔｉｏｎ）」，アンナルズ　オブ　インスティテュート　オブ　スタティスティカル
マセマティックス（Ａｎｎａｌｓ　ｏｆ　Ｉｎｓｔｉｔｕｔｅ　ｏｆ　Ｓｔａｔ．　Ｍａｔｈ．），１９６９年
【非特許文献２】
得丸・竹安：”離散時間線形モデルのあてはめによるスペクトル密度の推定”，計測自動制御学会論文集　１９７７年　１３−２　ｐ．４４−４９
【００１８】
【発明が解決しようとする課題】
（７）式に示す如きＡＲモデルに従う時系列｛ｙｔ　｝は、ｙｔ　の自己相関関数を、
Ｒｋ　＝Ｅ［ｙｔ　ｙｔ＋ｋ　］
とすると、以下のユール・ウォーカー方程式を満たす。
【００１９】
【数１７】
2004118355.tif 000009
【００２０】
時系列｛ｙｔ　｝が定常確率過程である場合は、ＡＲモデルのパラメータ｛αｍ　｝の推定値は、（９）式のユール・ウォーカー方程式の解が最良推定であることが知られている。このため、従来の方法では、（９）式を用いて前記パラメータを推定し、パワースペクトル密度の推定を行っていた。しかし、株価の変動または機械の振動などのデータからなる一般の時系列は、定常確率過程という仮定が成り立たず、（９）式を用いてパラメータを推定する方法では最適な推定ができないという問題があった。
【００２１】
本発明は、斯かる事情に鑑みてなされたものであって、その目的とするところは、ユール・ウォーカー方程式に加えて拡張ユール・ウォーカー方程式を用い、より多くのデータを利用してパワースペクトル密度の推定を行うことにより、より適したパワースペクトル密度の推定を行うことができるパワースペクトル密度推定方法、パワースペクトル密度推定装置、コンピュータを該パワースペクトル密度推定装置として実現するためのコンピュータプログラム、及び記録媒体を提供することにある。
【００２２】
【課題を解決するための手段】
第１発明に係るパワースペクトル密度推定方法は、記憶部及び演算部を備えたコンピュータを用いて、複数のデータを含む時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定する方法において、時系列をなす複数のデータを記憶部に記憶し、拡張ユール・ウォーカー方程式を用いて導かれる式により、記憶している時系列からＡＲモデルのパラメータの推定値を演算部にて計算し、計算した前記パラメータの推定値を用いて、ユール・ウォーカー方程式及び拡張ユール・ウォーカー方程式を用いて導かれる式により、前記時系列のパワースペクトル密度の推定値を演算部にて計算することを特徴とする。
【００２３】
第２発明に係るパワースペクトル密度推定方法は、解析対象から得られる時系列のデータを用いて前記解析対象の状態を解析するために、記憶部及び演算部を備えたコンピュータを用いて、前記時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定する方法において、解析対象から得られる時系列のデータ、及び該データが得られるときのサンプリングの時間間隔を記憶部に記憶し、拡張ユール・ウォーカー方程式を用いて導かれる式に、記憶している時系列のデータを代入して、ＡＲモデルのパラメータの推定値を演算部にて計算し、ユール・ウォーカー方程式及び拡張ユール・ウォーカー方程式を用いて導かれる式に、計算した前記パラメータの推定値、前記時系列のデータ、及び記憶している前記時間間隔を代入して、前記時系列のパワースペクトル密度の推定値を演算部にて計算することを特徴とする。
【００２４】
第３発明に係るパワースペクトル密度推定方法は、記憶部及び演算部を備えたコンピュータを用いて、複数のデータを含む時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定する方法において、時系列をなす複数のデータ｛ｙｔ　｝、及び該データ間の時間間隔Δを記憶部に記憶し、記憶している時系列の自己相関関数の推定値｛Ｒｋ　：ｋ＝０，１，…，Ｋ｝を、Ｒｋ　＝Ｅ［ｙｔ　ｙｔ＋ｋ　］と演算部にて計算し、前記時系列にＡＲモデルを適用したときのパラメータの推定値｛αｍ　：ｍ＝１，２，…，Ｍ｝を、
【００２５】
【数１８】
2004118355.tif 000010
【００２６】
と演算部にて計算し、記憶している前記時間間隔Δ及び計算したパラメータの推定値｛αｍ　｝を用いて、前記時系列のパワースペクトル密度の角周波数ωについての近似式Ｓ（ω）を、
【００２７】
【数１９】
2004118355.tif 000011
【００２８】
の式を用いて演算部にて計算することを特徴とする。
【００２９】
第４発明に係るパワースペクトル密度推定方法は、記憶部及び演算部を備えたコンピュータを用いて、複数のデータを含む時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定する方法において、時系列をなす複数のデータ｛ｙｔ　｝、及び該データ間の時間間隔Δを記憶部に記憶し、記憶している時系列の自己相関関数の推定値｛Ｒｋ　：ｋ＝０，１，…，Ｋ｝を、Ｒｋ　＝Ｅ［ｙｔ　ｙｔ＋ｋ　］と演算部にて計算し、前記時系列にＡＲモデルを適用したときのパラメータの推定値｛αｍ　：ｍ＝１，２，…，Ｍ｝を、
【００３０】
【数２０】
2004118355.tif 000012
【００３１】
の式を用いて演算部にて計算し、記憶している前記時間間隔Δ及び計算したパラメータの推定値｛αｍ　｝を用いて、前記時系列のパワースペクトル密度の角周波数ωについての近似式Ｓ（ω）を、
【００３２】
【数２１】
2004118355.tif 000013
【００３３】
と演算部にて計算することを特徴とする。
【００３４】
第５発明に係るパワースペクトル密度推定装置は、複数のデータを含む時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定する装置において、時系列をなす複数のデータ｛ｙｔ　｝、及び該データ間の時間間隔Δを記憶する手段と、記憶している時系列の自己相関関数の推定値｛Ｒｋ　：ｋ＝０，１，…，Ｋ｝を、Ｒｋ　＝Ｅ［ｙｔ　ｙｔ＋ｋ　］と計算する手段と、前記時系列にＡＲモデルを適用したときのパラメータの推定値｛αｍ　：ｍ＝１，２，…，Ｍ｝を、
【００３５】
【数２２】
2004118355.tif 000014
【００３６】
の式を用いて計算する手段と、記憶している前記時間間隔Δ及び計算したパラメータの推定値｛αｍ　｝を用いて、前記時系列のパワースペクトル密度の角周波数ωについての近似式Ｓ（ω）を、
【００３７】
【数２３】
2004118355.tif 000015
【００３８】
と計算する手段とを備えることを特徴とする。
【００３９】
第６発明に係るコンピュータプログラムは、コンピュータに、複数のデータを含む時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定させるコンピュータプログラムにおいて、コンピュータに、記憶している時系列をなす複数のデータ｛ｙｔ　｝の自己相関関数の推定値｛Ｒｋ　：ｋ＝０，１，…，Ｋ｝を、Ｒｋ　＝Ｅ［ｙｔ　ｙｔ＋ｋ　］と計算させる手順と、コンピュータに、前記時系列にＡＲモデルを適用したときのパラメータの推定値｛αｍ　：ｍ＝１，２，…，Ｍ｝を、
【００４０】
【数２４】
2004118355.tif 000016
【００４１】
の式を用いて計算させる手順と、コンピュータに、記憶している前記データ間の時間間隔Δ、及び計算したパラメータの推定値｛αｍ　｝を用いて、前記時系列のパワースペクトル密度の角周波数ωについての近似式Ｓ（ω）を、
【００４２】
【数２５】
2004118355.tif 000017
【００４３】
と計算させる手順とを含むことを特徴とする。
【００４４】
第７発明に係るコンピュータでの読み取りが可能な記録媒体は、コンピュータに、複数のデータを含む時系列にＡＲモデルを適用して前記時系列のパワースペクトル密度を推定させるコンピュータプログラムを記録してあるコンピュータでの読み取りが可能な記録媒体において、コンピュータに、記憶している時系列をなす複数のデータ｛ｙｔ　｝の自己相関関数の推定値｛Ｒｋ　：ｋ＝０，１，…，Ｋ｝を、Ｒｋ　＝Ｅ［ｙｔ　ｙｔ＋ｋ　］と計算させる手順と、コンピュータに、前記時系列にＡＲモデルを適用したときのパラメータの推定値｛αｍ　：ｍ＝１，２，…，Ｍ｝を、
【００４５】
【数２６】
2004118355.tif 000018
【００４６】
の式を用いて計算させる手順と、コンピュータに、記憶している前記データ間の時間間隔Δ、及び計算したパラメータの推定値｛αｍ　｝を用いて、前記時系列のパワースペクトル密度の角周波数ωについての近似式Ｓ（ω）を、
【００４７】
【数２７】
2004118355.tif 000019
【００４８】
と計算させる手順とを含むコンピュータプログラムを記録してあることを特徴とする。
【００４９】
本発明者は、時系列にＡＲモデルを適用して該時系列のパワースペクトル密度を推定する際に、ユール・ウォーカー方程式に加えて拡張ユール・ウォーカー方程式を用いてパワースペクトル密度の近似式を導出する事で、より良いパワースペクトル密度の推定を試みた。（７）式に示す如きＡＲモデルに従う時系列｛ｙｔ　｝は、以下の拡張ユール・ウォーカー方程式を満たす。
【００５０】
【数２８】
2004118355.tif 000020
【００５１】
Ｒｎ　＝Ｒ−ｎであるため、α０　＝１とすることにより、（９）式および（１０）式は、
【００５２】
【数２９】
2004118355.tif 000021
【００５３】
と表すことができる。（９）式に含まれるＭ個の式の夫々について、１番目の式から順次ｅ−ｉωｎΔ　（ｎ＝１，２，…，Ｍ）を掛け、（１０）式に含まれる複数の式の夫々について、第１番目の式から順次ｅ−ｉωｎΔ　（ｎ＝Ｍ＋１，Ｍ＋２，…）を掛け、更にｅ−ｉωｎΔ　を掛けた夫々の式の全てについて和をとると、（１１）式より、
【００５４】
【数３０】
2004118355.tif 000022
【００５５】
の式を得ることができる。ここで、
【００５６】
【数３１】
2004118355.tif 000023
【００５７】
とおくと、（１２）式に（１３）式を代入して変形することにより、
【００５８】
【数３２】
2004118355.tif 000024
【００５９】
となる。（３）式と（１３）式とを比較することにより、時系列｛ｙｔ　｝のパワースペクトルＳ（ω）の近似式は、
Ｓ（ω）＝Δ｛Ｑ（ω）＋Ｑ（ω）＊　＋Ｒ０　｝　　　　…（１５）
となる。ＡＲモデルのパラメータが得られた場合は、（１４）式および（１５）式を用いることにより、パワースペクトル密度の推定値を計算できる。
【００６０】
次に、ＡＲモデルのパラメータ｛αｍ　：ｍ＝１，２，…，Ｍ｝の推定を試みる。時系列｛ｙｔ　｝の自己相関関数の統計的誤差を含んだ推定値｛Ｒｋ　：ｋ＝０，１，…，Ｋ｝が得られたときに、拡張ユール・ウォーカー方程式を利用して｛αｍ　｝の推定値を求める。拡張ユール・ウォーカー方程式の右辺は（１０）式の如く理論的には０となるが、｛Ｒｋ　｝が統計的誤差を含んでいるために誤差εが生じると考えて、（１０）式を、
ＲＭ＋ｎ　＋α１　ＲＭ＋ｎ−１　＋…＋αＭ　Ｒｎ　＝εｎ　　　（ｎ＝１，２，…，Ｌ）…（１６）
とおく。但し、Ｌ＋Ｍ≦Ｋとする。（１６）式のεｎ　の２乗和が最小になるように｛αｍ　｝を定める。ここで、行列α（Ｍ）及び行列Ａｎ　を、
α（Ｍ）＝［αＭ　，αＭ−１　，…，α１　］Ｔ
Ａｎ　＝［Ｒｎ　，Ｒｎ＋１　，…，ＲＭ＋ｎ−１　］
と定義すると、（１６）式は、
ＲＭ＋ｎ　＋Ａｎ　α（Ｍ）＝εｎ　　　　　　　　　…（１７）
となる。（１７）式の２乗和を取ると、
【００６１】
【数３３】
2004118355.tif 000025
【００６２】
となる。｛αｍ　｝を定めるためには、（１８）式を最小にするα（Ｍ）を求めればよい。（１８）式をα（Ｍ）で偏微分して０に等しいとおくことにより、
【００６３】
【数３４】
2004118355.tif 000026
【００６４】
となる。ここで、行列ｒｎ　、行列Ｒ（Ｍ）、行列Ｆ（Ｍ）及び行列ｆ（ｍ）を、
ｒｎ　＝［Ｒｎ　，Ｒｎ＋１　，…，Ｒｎ−Ｌ　］Ｔ
Ｒ（Ｍ）＝［ｒ０　，ｒ１　，…，ｒＭ−１　］
Ｆ（Ｍ）＝Ｒ（Ｍ）Ｔ　Ｒ（Ｍ）
ｆ（Ｍ）＝Ｒ（Ｍ）Ｔ　ｒＭ
と定義すると、（１９）式は、
Ｒ（Ｍ）Ｔ　Ｒ（Ｍ）α（Ｍ）＝−Ｒ（Ｍ）Ｔ　ｒＭ
となり、
α（Ｍ）＝−Ｆ（Ｍ）−１ｆ（Ｍ）　　…（２０）
が得られる。ここで、Ｆ（Ｍ）−１については、
ｇ（Ｍ）＝ｒＭ　Ｔ　ｒＭ
を定義すると、
【００６５】
【数３５】
2004118355.tif 000027
【００６６】
の関係がある。（２０）式および（２１）式を用いることで、Ｍ＝１からＭを逐次増加させて、逆行列の計算を行うことなくα（Ｍ）を計算することができる。
【００６７】
第１発明においては、時系列のパワースペクトル密度を推定する際に、前記時系列にＡＲモデルを適用したときのパラメータの推定値を、拡張ユール・ウォーカー方程式を利用して計算し、ユール・ウォーカー方程式及び拡張ユール・ウォーカー方程式を利用して、前記時系列と前記パラメータの推定値とから前記パワースペクトル密度の近似式を計算することにより、一般の非定常確率過程の時系列についても、ユール・ウォーカー方程式のみを用いてパワースペクトル密度を推定していた従来の方法に比べてより良いパワースペクトル密度の推定が可能となる。
【００６８】
第２発明においては、スペクトル解析の解析対象から得られる時系列のデータからパワースペクトル密度を推定する際に、ユール・ウォーカー方程式及び拡張ユール・ウォーカー方程式を利用して、前記解析対象の状態を示す時系列のデータ、及び該データのサンプリング間隔から、前記パワースペクトル密度の近似式を計算することにより、従来の方法に比べてより良いパワースペクトル密度の推定が可能となる。
【００６９】
第３発明においては、時系列のパワースペクトル密度を推定する際に、前記時系列にＡＲモデルを適用したときのパラメータの推定値を、拡張ユール・ウォーカー方程式を利用して導かれた（２０）式を用いて計算し、ユール・ウォーカー方程式及び拡張ユール・ウォーカー方程式を利用して導かれた（１５）式を用いて前記パワースペクトル密度の近似式を計算することにより、一般の時系列について、従来の方法に比べてより良いパワースペクトル密度の近似式が得られる。
【００７０】
第４、第５、第６及び第７発明においては、時系列のパワースペクトル密度を推定する際に、前記時系列にＡＲモデルを適用したときのパラメータの推定値を、拡張ユール・ウォーカー方程式を利用して導かれた（２０）式及び（２１）式を用いて計算し、ユール・ウォーカー方程式及び拡張ユール・ウォーカー方程式を利用して導かれた（１５）式を用いて前記パワースペクトル密度の近似式を計算することにより、一般の時系列について、従来の方法に比べてより良いパワースペクトル密度の近似式が容易に得られる。
【００７１】
【発明の実施の形態】
以下本発明をその実施の形態を示す図面に基づき具体的に説明する。
図１は、本発明のパワースペクトル密度推定装置の構成を示すブロック図である。パワースペクトル密度推定装置１は、コンピュータを用いて構成されており、演算を行うＣＰＵ（演算部）１１と、演算に伴って発生する一時的な情報を記憶するＲＡＭ（記憶部）１２と、ＣＤ−ＲＯＭドライブ等の外部記憶装置１３と、ハードディスク等の内部記憶装置１４とを備えており、ＣＤ−ＲＯＭ等の本発明に係る記録媒体２から本発明に係るコンピュータプログラム２０を外部記憶装置１３にて読み取り、読み取ったコンピュータプログラム２０を内部記憶装置１４に記憶し、ＲＡＭ１２にコンピュータプログラム２０をロードし、ＣＰＵ１１は、ロードしたコンピュータプログラム２０に基づいてパワースペクトル密度推定装置１に必要な処理を実行する。パワースペクトル密度推定装置１は、キーボード又はマウス等の入力装置１５と、液晶ディスプレイ又はＣＲＴディスプレイ等の出力装置１６とを備えており、データの入力を初めとするオペレータからの操作を受け付ける構成となっている。
【００７２】
なお、パワースペクトル密度推定装置１は、図示しない通信インタフェースを備え、通信インタフェースに接続されている図示しない外部のサーバ装置から本発明に係るコンピュータプログラム２０をダウンロードし、ＣＰＵ１１にて処理を実行する形態であってもよい。
【００７３】
次に、フローチャートを用いて本発明のパワースペクトル密度推定方法の手順を説明する。図２は、本発明のパワースペクトル密度推定装置１が行う処理の手順を示すフローチャートである。パワースペクトル密度推定装置１は、オペレータの操作により入力装置１５に入力されたデータを読み込む、又は内部記憶装置１４に記憶してあるデータを読み込む等して、Ｎ個のデータを含む時系列｛ｙｔ　：ｔ＝１，２，…｝、及び該時系列におけるデータの時間間隔ΔをＲＡＭ１２に読み込む（Ｓ１）。時系列のデータは、株価もしくは販売量などのデータ、又は振動もしくは温度などの機械の状態を示すデータ等、定常確率過程が保証されない時系列のデータであり、時間間隔は、前記データのサンプリング間隔である。ＣＰＵ１１は、コンピュータプログラム２０をＲＡＭ１２へロードし、ロードしたコンピュータプログラム２０に従って、読み込んだ時系列の自己相関関数の推定値｛Ｒｋ　｝を該時系列から計算する（Ｓ２）。ＣＰＵ１１は、次に、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、Ｍ＝１とし（Ｓ３）、（２０）式を用いてα（Ｍ）を計算し（Ｓ４）、計算したα（Ｍ）をＲＡＭ１２に記憶する（Ｓ６）。ＣＰＵ１１は、次に、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、Ｍの値が所定の値に達しているか否かを判定し（Ｓ６）、Ｍの値が小さくて所定値に達していない場合は（Ｓ６：ＮＯ）、（２１）式を用いてＦ（Ｍ＋１）−１を計算し（Ｓ７）、Ｍの値を１増加させ（Ｓ８）、処理をステップＳ４へ戻してα（Ｍ）を計算する。
【００７４】
ステップＳ６にてＭが所定値に達していた場合は（Ｓ６：ＹＥＳ）、ＣＰＵ１１は、次に、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、夫々のＭの値について計算したα（Ｍ）の値を利用した非特許文献１に開示されたＦＰＥ法または非特許文献２に開示されたＦＯＥ法を用いて、時系列｛ｙｔ　｝にＡＲモデルを適用したときの最適なＡＲモデルの次数Ｍの値を決定する（Ｓ９）。ＣＰＵ１１は、次に、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、ＲＡＭ１２に記憶してある夫々のＭの値に対応したα（Ｍ）から、決定したＭの値に対応したα（Ｍ）を抜き出して、ＡＲモデルのパラメータの推定値である｛αｍ　｝を決定し（Ｓ１０）、決定した｛αｍ　｝の値を利用して、（１４）式よりＱ（ω）を計算し（Ｓ１１）、計算したＱ（ω）から、（１５）式を用いてパワースペクトル密度の近似式Ｓ（ω）を計算し（Ｓ１２）、処理を終了する。
【００７５】
なお、前述の処理では、（２１）式を用いてＦ（Ｍ＋１）−１を計算する処理を用いているが、これに限る物ではなく、逆行列の計算をしてＦ（Ｍ）−１を直接に求め、求めたＦ（Ｍ）−１を用いて（２０）式によりα（Ｍ）を計算する処理を用いてもよい。また、前述の処理では、Ｍが所定値に達するまでの全てのＭについてα（Ｍ）を計算し、ＦＰＥ法またはＦＯＥ法を用いて最適なＭを決定する処理を用いているが、α（Ｍ）の値を求めるためのＭの値を所定のモデルに従って予め定めておき、定めてあるＭの値についてα（Ｍ）の値を求める処理を用いても良い。
【００７６】
次に、本発明と従来技術とで推定されたパワースペクトル密度を比較する。パワースペクトル密度が既知である時系列｛ｙｔ　｝を作成し、作成した｛ｙｔ　｝について本発明および従来技術を用いてパワースペクトル密度の推定を行った。図３は、本発明により推定されたパワースペクトル密度と従来技術により推定されたパワースペクトル密度とを示す特性図である。図の縦軸はパワースペクトル密度の値をログスケールで示し、横軸は角周波数ωを２πで除した値をログスケールで示している。図中には、実際のパワースペクトル密度を実線で示し、本発明により推定されたパワースペクトル密度を一点鎖線で示し、従来技術により推定されたパワースペクトル密度を破線で示している。本発明では、従来技術に比べて、パワースペクトル密度の再現性に優れている。これにより、ユール・ウォーカー方程式に加えて拡張ユール・ウォーカー方程式を用いることで、ユール・ウォーカー方程式のみを用いる従来の方法に比べて、より多くの情報に基づいて最適値を求めるため、ユール・ウォーカー方程式の解が最良の推定値を与える保証がない一般の時系列について、より正確なパワースペクトル密度の推定が可能となっていることがわかる。また、図３により、特に低周波領域のパワースペクトル密度をより正しく推定することが明らかであり、これは、本発明が、時系列をなすデータの中で従来技術よりもより離れたデータ間の相関関係をも利用してパワースペクトル密度を推定しているからである。
【００７７】
以上詳述した如く、本発明においては、時系列をなす複数のデータについて、従来技術に比べてより良いパワースペクトル密度の推定を行うことができる。また、従来技術に比べて、より低周波、即ちより長周期のパワースペクトル密度をより正しく推定することができる。これにより、本発明を振動などの観測データの分析に用いる場合は、低周波の振動モードの検出など低周波領域での現象を適切に解析することが可能となり、また、本発明を数値の変動予測に利用する場合は、長周期での数値の変動を予測することが可能となる。
【００７８】
【発明の効果】
第１発明においては、時系列のパワースペクトル密度を推定する際に、ユール・ウォーカー方程式のみを用いてパワースペクトル密度を推定していた従来の方法に比べ、一般の時系列についてより良いパワースペクトル密度の推定が可能となる。
【００７９】
第２発明においては、スペクトル解析の解析対象から得られる時系列のデータからパワースペクトル密度を推定する際に、ユール・ウォーカー方程式及び拡張ユール・ウォーカー方程式を利用して、前記解析対象の状態を示す時系列のデータ、及び該データのサンプリング間隔から、前記パワースペクトル密度の近似式を計算することにより、従来の方法に比べてより良いパワースペクトル密度の推定が可能となる。また、従来の方法に比べてより低周波のパワースペクトル密度をより正しく推定することができ、低周波領域での現象の解析または長周期での数値の変動予測などが可能となる。
【００８０】
第３発明においては、時系列のパワースペクトル密度を推定する際に、前記時系列にＡＲモデルを適用したときのパラメータの推定値を、拡張ユール・ウォーカー方程式を利用して導かれた（２０）式を用いて計算し、ユール・ウォーカー方程式及び拡張ユール・ウォーカー方程式を利用して導かれた（１５）式を用いて前記パワースペクトル密度の近似式を計算することにより、一般の時系列について、従来の方法に比べてより良いパワースペクトル密度の近似式が得られる。
【００８１】
第４、第５、第６及び第７発明においては、時系列のパワースペクトル密度を推定する際に、前記時系列にＡＲモデルを適用したときのパラメータの推定値を、拡張ユール・ウォーカー方程式を利用して導かれた（２０）式及び（２１）式を用いて計算し、ユール・ウォーカー方程式及び拡張ユール・ウォーカー方程式を利用して導かれた（１５）式を用いて前記パワースペクトル密度の近似式を計算することにより、一般の時系列について、従来の方法に比べてより良いパワースペクトル密度の近似式が容易に得られる等、本発明は優れた効果を奏する。
【図面の簡単な説明】
【図１】本発明のパワースペクトル密度推定装置の構成を示すブロック図である。
【図２】本発明のパワースペクトル密度推定装置が行う処理の手順を示すフローチャートである。
【図３】本発明により推定されたパワースペクトル密度と従来技術により推定されたパワースペクトル密度とを示す特性図である。
【符号の説明】
１　パワースペクトル密度推定装置
１１　ＣＰＵ（演算部）
１２　ＲＡＭ（記憶部）
２　記録媒体
２０　コンピュータプログラム
｛ｙｔ　｝　時系列
｛αｍ　｝　ＡＲモデルのパラメータ
Ｓ（ω）　パワースペクトル密度
【図１】
2004118355.tif 000038
【図２】
2004118355.tif 000039
【図３】
2004118355.tif 000040