(57)【要約】【課題】パラメータの値の収束を容易にし、同時に計算時間の短縮を図ることができるパラメータ推定方法を提供する。【解決手段】３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式を用いて、まずＡＲパラメータであるベクトルａを計算する（ステップＳ２３）。さらに白色雑音｛ｅｎ  ｝を、乱数等を用いて生成し（ステップＳ２４）、生成した白色雑音と記憶部に記憶した時系列とから、ベクトルｂを計算する（ステップＳ２７）。演算部は、計算したベクトルａ及びｂから、新たな白色雑音｛ｅｎ  ｝を計算し（ステップＳ２９）、推定すべきパラメータであるベクトルｂが収束するまで計算を繰り返す。【選択図】        図２
2004302690.tif 000002
【特許請求の範囲】
【請求項１】
記憶部及び演算部を備えたコンピュータを用いて、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときのパラメータを推定する方法において、
Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を記憶部に記憶し、
平均値が０であるＮ個の値が時系列で並んだ白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を演算部にて生成し、
生成した白色雑音の分散σｅ２を演算部にて計算し、
前記ＡＲＭＡモデルにおけるＡＲ過程のｐ個のパラメータを行列要素とする行列ａ＝［ａ１，ａ２，…，ａｐ］Ｔを３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式を用いて前記演算部により計算し、
前記ＡＲＭＡモデルにおけるＭＡ過程のｐ個のパラメータを行列要素とする行列ｂ＝［ｂ１，ｂ２，…，ｂｐ］Ｔを、前記計算により求めた行列ａ及び数１を用いて、
【数１】
2004302690.tif 000025
と演算部にて計算し、
計算した行列ａ及び行列ｂに含まれる前記ＡＲＭＡモデルのパラメータ、並びに前記時系列の値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を演算部にて計算し、
計算した新たな白色雑音を用いて新たな行列ａ及び行列ｂを演算部にて計算し、
行列ａ及び行列ｂが収束するまで行列ａ及び行列ｂの計算を演算部にて繰り返す
ことを特徴とするパラメータ推定方法。
【請求項２】
行列ａ＝［ａ１，ａ２，…，ａｐ］Ｔは、数２に基づいて前記計算部により計算される
【数２】
2004302690.tif 000026
ことを特徴とする請求項１に記載のパラメータ推定方法。
【請求項３】
記憶部及び演算部を備えたコンピュータを用いて、複数のデータから構成される時系列から該時系列に加わる新たなデータを予測する方法において、
時系列を記憶部に記憶し、
請求項１または２に記載のパラメータ推定方法を用い、前記時系列をＡＲＭＡモデルに当てはめたときのパラメータを演算部にて計算し、
計算したパラメータ、及び記憶している前記時系列の値を前記ＡＲＭＡモデルの式に代入して、前記時系列に加わる新たなデータを演算部にて計算する
ことを特徴とするデータ予測方法。
【請求項４】
Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときのパラメータを推定する装置において、
Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を記憶する手段と、
平均値が０であるＮ個の値が時系列で並んだ白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成する手段と、
生成した白色雑音の分散σｅ２を計算する手段と、
前記ＡＲＭＡモデルにおけるＡＲ過程のｐ個のパラメータを行列要素とする行列ａ＝［ａ１，ａ２，…，ａｐ］Ｔを３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式を用いて前記演算部により計算する手段と、
前記ＡＲＭＡモデルにおけるＭＡ過程のｐ個のパラメータを行列要素とする行列ｂ＝［ｂ１，ｂ２，…，ｂｐ］Ｔを、前記計算により求めた行列ａ及び数３を用いて、
【数３】
2004302690.tif 000027
と演算部にて計算する手段と、
計算した行列ａ及び行列ｂに含まれる前記ＡＲＭＡモデルのパラメータ、並びに前記時系列の値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算する手段と、
計算した新たな白色雑音を用いて新たな行列ａ及び行列ｂを計算する手段と、
行列ａ及び行列ｂが収束するまで行列ａ及び行列ｂの計算を繰り返す手段と
を備えることを特徴とするパラメータ推定装置。
【請求項５】
複数のデータから構成される時系列から該時系列に加わる新たなデータを予測する装置において、
時系列を記憶する手段と、
請求項４に記載のパラメータ推定装置と、
該パラメータ推定装置を用いて、前記時系列から、前記時系列をＡＲＭＡモデルに当てはめたときのパラメータを計算する手段と、
計算したパラメータ、及び記憶している前記時系列の値を前記ＡＲＭＡモデルの式に代入して、前記時系列に加わる新たなデータを計算する手段と
を備えることを特徴とするデータ予測装置。
【請求項６】
コンピュータに、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときのパラメータを推定させるコンピュータプログラムにおいて、
コンピュータに、平均値が０であるＮ個の値が時系列で並んだ白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成させる手順と、
コンピュータに、生成した白色雑音の分散σｅ２を計算させる手順と、
コンピュータに、記憶している、Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を用いて、前記ＡＲＭＡモデルにおけるＡＲ過程のｐ個のパラメータを行列要素とする行列ａ＝［ａ１，ａ２，…，ａｐ］Ｔを３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式により計算させる手順と、
コンピュータに、記憶している、Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を用いて、前記ＡＲＭＡモデルにおけるＭＡ過程のｐ個のパラメータを行列要素とする行列ｂ＝［ｂ１，ｂ２，…，ｂｐ］Ｔを、前記計算により求めた行列ａ及び数４により、
【数４】
2004302690.tif 000028
と計算させる手順と、
コンピュータに、計算した行列ａ及び行列ｂに含まれる前記ＡＲＭＡモデルのパラメータ、並びに前記時系列の値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算させる手順と、
コンピュータに、計算した新たな白色雑音を用いて新たな行列ａ及び行列ｂを計算させる手順と、
コンピュータに、行列ａ及び行列ｂが収束するまで行列ａ及び行列ｂの計算を繰り返させる手順と
を含むことを特徴とするコンピュータプログラム。
【請求項７】
コンピュータに、複数のデータから構成される時系列から該時系列に加わる新たなデータを予測させるコンピュータプログラムにおいて、
コンピュータに、請求項６に記載のコンピュータプログラムを用いて、記憶している時系列から、該時系列をＡＲＭＡモデルに当てはめたときのパラメータを計算させる手順と、
コンピュータに、計算したパラメータ、及び記憶している前記時系列の値を前記ＡＲＭＡモデルの式に代入して、前記時系列に加わる新たなデータを計算させる手順と
を含むことを特徴とするコンピュータプログラム。
【請求項８】
コンピュータに、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときのパラメータを推定させるコンピュータプログラムを記録してあるコンピュータでの読み取りが可能な記録媒体において、
コンピュータに、平均値が０であるＮ個の値が時系列で並んだ白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成させる手順と、
コンピュータに、生成した白色雑音の分散σｅ２を計算させる手順と、
コンピュータに、記憶している、Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を用いて、前記ＡＲＭＡモデルにおけるＡＲ過程のｐ個のパラメータを行列要素とする行列ａ＝［ａ１，ａ２，…，ａｐ］Ｔを３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式により計算させる手順と、
コンピュータに、記憶している、Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を用いて、前記ＡＲＭＡモデルにおけるＭＡ過程のｐ個のパラメータを行列要素とする行列ｂ＝［ｂ１，ｂ２，…，ｂｐ］Ｔを、前記計算により求めた行列ａ及び数５により、
【数５】
2004302690.tif 000029
と計算させる手順と、
コンピュータに、計算した行列ａ及び行列ｂに含まれる前記ＡＲＭＡモデルのパラメータ、並びに前記時系列の値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算させる手順と、
コンピュータに、計算した新たな白色雑音を用いて新たな行列ａ及び行列ｂを計算させる手順と、
コンピュータに、行列ａ及び行列ｂが収束するまで行列ａ及び行列ｂの計算を繰り返させる手順と
を含むコンピュータプログラムを記録してあることを特徴とするコンピュータでの読み取りが可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、データの時系列からパラメータを推定する方法、推定されたパラメータを用いて新たなデータを予測する方法、パラメータ推定装置、データ予測方法、コンピュータをパラメータ推定装置またはデータ予測装置として実現するためのコンピュータプログラム、及び記録媒体に関する。
【０００２】
【従来の技術】
蓄積されたデータの時系列を解析し、新たなデータの値を予想する方法は、従来、様々な方法が提案されている。株価または販売量などのデータから構成される一般の経済的な時系列は、白色ノイズによって励起された系の出力とみることができる。データの時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ，…｝を、定常エルゴード的正規過程ｘ（ｔ）の標本時系列であるとしたときには、（ｐ，ｑ）次のＡＲＭＡモデルでは、｛ｘｎ  ｝の要素間の関係は、以下の式で表される
【０００３】
【数６】
2004302690.tif 000003
【０００４】
ここで、｛ｅｎ  ｝は、平均値０及び分散σｅ２の正規性雑音であり、｛ａｉ  ：ｉ＝１，２，…，ｐ｝及び｛ｂｊ  ：ｊ＝１，２，…，ｑ｝は、
Ａ（Ｚ−１）＝１＋ａ１  Ｚ−１＋…＋ａｐ  Ｚ−ｐ
Ｂ（Ｚ−１）＝１＋ｂ１  Ｚ−１＋…＋ｂｑ  Ｚ−ｑ
とおいたときに、Ａ（Ｚ−１）及びＢ（Ｚ−１）が規約であり、更に定常条件、可逆条件および強正実条件を満たすものとする。
【０００５】
標本の時系列｛ｘｎ  ｝から（１）式のパラメータ｛ａｉ  ｝及び｛ｂｊ  ｝を推定するための方法として、非特許文献１に開示された疑似線形回帰法がある。簡単にするためにｐ＝ｑとし、ベクトルＺｎ  及びベクトルθを次の様に定義する。
Ｚｎ  ＝［−ｘｎ−１，…，−ｘｎ−ｐ，ｅｎ−１，…，ｅｎ−ｐ］Ｔ
θ＝［ａ１，ａ２，…，ａｐ，ｂ１，ｂ２，…，ｂｐ］Ｔ
但し、ｂｐ  ＝０
（１）式をｘｎ  についての式に書き直したときには、（１）式は、
ｘｎ  ＝θＴ  Ｚｎ  ＋ｅｎ      …（２）
と表すことができる。Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝が得られたとき、ベクトルθを時系列から最小二乗法により求めるためには、
【０００６】
【数７】
2004302690.tif 000004
【０００７】
の右辺を最小にするθを求めればよい。従って、（３）式の右辺をθで偏微分して０に等しいとすることにより、時系列から推定したベクトルθは、次の式で与えられる。
【０００８】
【数８】
2004302690.tif 000005
【０００９】
ラグがｋである｛ｘｎ  ｝の自己相関関数をＲｋ  とすると、
Ｒｋ  ＝Ｅ［ｘｎ  ｘｎ＋ｋ  ］
である。ここで、Ｅ［ｙｎ  ］は、データ列｛ｙｎ  ：ｎ＝１，２，…｝の期待値である。また、
Ｒ−ｋ＝Ｒｋ
である。ラグがｋである｛ｅｎ  ｝の自己相関関数をＳｋとすると、同様に、
Ｓｋ＝Ｅ［ｅｎ  ｅｎ＋ｋ  ］，Ｓ−ｋ＝Ｓｋ
である。また、｛ｘｎ  ｝と｛ｅｎ  ｝との相互相関関数は、
Ｔｋ  ＝Ｅ［ｅｎ  ｘｎ＋ｋ  ］
Ｔ−ｋ＝Ｅ［ｅｎ  ｘｎ−ｋ  ］
とおくと、
Ｅ［ｘｎ  ｅｎ＋ｋ  ］＝Ｅ［ｅｎ  ｘｎ−ｋ  ］＝Ｔ−ｋ
Ｅ［ｘｎ  ｅｎ−ｋ  ］＝Ｅ［ｅｎ  ｘｎ＋ｋ  ］＝Ｔｋ
となる。
【００１０】
（４）式の右辺は、｛ｘｎ  ｝の自己相関関数、｛ｅｎ  ｝の自己相関関数、及び｛ｘｎ  ｝と｛ｅｎ  ｝との相互相関関数で表すことができ、（４）式は、
【００１１】
【数９】
2004302690.tif 000006
【００１２】
となる。ベクトルａ及びｂを、
ａ＝［ａ１，ａ２，…，ａｐ］Ｔ
ｂ＝［ｂ１，ｂ２，…，ｂｐ］Ｔ
と定義し、（５）式に含まれる行列を、以下の（６）式と定義したときには、（５）式は（７）式と表すことができる。
【００１３】
【数１０】
2004302690.tif 000007
【００１４】
ベクトルａ，ｂ，ｒ，ｔはｐ行１列の行列であり、行列Ｒ，Ｔ，Ｓはｐ行ｐ列の行列であるため、（７）式はベクトルａ及びｂについての連立方程式の形に書き直すことができ、
Ｒ・ａ−τ・ｂ＝−ｒ      …（８）
−τＴ  ・ａ＋Ｓ・ｂ＝ｔ    …（９）
となる。（９）式より、
ｂ＝Ｓ−１（ｔ＋τＴ  ・ａ）        …（１０）
であり、（１０）式を（８）式へ代入することにより、
ａ＝（Ｒ−τ・Ｓ−１・τＴ  ）−１（τ・Ｓ−１・ｔ−ｒ）    …（１１）
が求められる。（１１）式よりベクトルａが得られ、得られたベクトルａを（１０）式に代入することにより、ベクトルｂを計算することができる。
【００１５】
疑似線形回帰法では、白色雑音｛ｅｎ  ｝を乱数などを用いて生成し、生成した｛ｅｎ  ｝と時系列｛ｘｎ  ｝とから、（１１）式および（１０）式を用いて、（１）式のパラメータ｛ａｉ  ｝及び｛ｂｊ  ｝を計算し、計算した｛ａｉ  ｝及び｛ｂｊ  ｝から、（１）式を用いて白色雑音｛ｅｎ  ｝を計算し、パラメータ｛ａｉ  ｝及び｛ｂｊ  ｝が収束するまで計算を繰り返し、パラメータ｛ａｉ  ｝及び｛ｂｊ  ｝を求める。
【００１６】
また、ＡＲモデルのパラメータはＹｕｌｅ−Ｗａｌｋｅｒ方程式により求まることが公知であるところ、高次統計量キュムラントに対しても同様のＹｕｌｅ−Ｗａｌｋｅｒ方程式が得られることが知られている。以下に３次キュムラントについて考える。ＡＲ（ｐ）モデルは（１２）式で記述される。
【００１７】
【数１１】
2004302690.tif 000008
【００１８】
ここで、（１２）式のモデルは安定と仮定とする。いま、（１３）式と略記する。
【００１９】
【数１２】
2004302690.tif 000009
【００２０】
（１２）式の両辺にｙ（ｔ−ｍ）２を乗じて平均をとった場合、入力過程｛ｕ（ｔ）｝が３次白色であるので、（１４）式を得ることができる。
【００２１】
【数１３】
2004302690.tif 000010
【００２２】
（１４）式において、ｍを１からｐまで動かす場合（１５）式で示す行列方程式を得ることができる。
【００２３】
【数１４】
2004302690.tif 000011
【００２４】
（１５）式は３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式といわれる。（１５）式を用いることにより、パラメータ｛ａｉ  ｝を容易に求めることができる（例えば、非特許文献２、非特許文献３）。
【００２５】
【非特許文献１】
片山徹「システム同定入門」朝倉書店  １９９４
【非特許文献２】
井上雄二郎  「高次統計量キュムラントによる信号処理−Ｉ」  システム／制御／情報，Ｖｏｌ．３６，  Ｎｏ．２，  ｐｐ．９０〜９９，  １９９２
【非特許文献３】
井上雄二郎  「高次統計量キュムラントによる信号処理−ＩＩ」  システム／制御／情報，Ｖｏｌ．３６，  Ｎｏ．５，  ｐｐ．２９４〜３０６，  １９９２
【００２６】
【発明が解決しようとする課題】
しかしながら、疑似線形回帰法は、問題の性質によっては、パラメータ｛ａｉ  ｝及び｛ｂｊ  ｝が収束せず、パラメータの推定ができない場合があるという問題がある。また、ｐ行ｐ列の行列計算を行うため、計算時間が長いという問題があった。
【００２７】
本発明は斯かる事情に鑑みてなされたものであり、その目的とするところは、３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式を用いてＡＲ項のパラメータを推定し、ＭＡ項のパラメータについては計算に用いる数値のうち理論的に０となるはずの数値を０とおいて計算を行うことにより、パラメータの値の収束を容易にし、同時に計算時間の短縮を図ることができるパラメータ推定方法、パラメータ推定装置、コンピュータを該パラメータ推定装置として実現するためのコンピュータプログラム、及び記録媒体を提供することにある。
【００２８】
また、本発明の他の目的とするところは、前記パラメータ推定方法により推定したパラメータを用いてデータの推定を行うことにより、従来の手法ではデータの予測ができなかった時系列についてもデータの予測が可能となるデータ予測方法、データ予測装置、コンピュータを該データ予測装置として実現するためのコンピュータプログラム、及び記録媒体を提供することにある。
【００２９】
【課題を解決するための手段】
本発明に係るパラメータ推定方法は、記憶部及び演算部を備えたコンピュータを用いて、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときのパラメータを推定する方法において、Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を記憶部に記憶し、平均値が０であるＮ個の値が時系列で並んだ白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を演算部にて生成し、生成した白色雑音の分散σｅ２を演算部にて計算し、前記ＡＲＭＡモデルにおけるＡＲ過程のｐ個のパラメータを行列要素とする行列ａ＝［ａ１，ａ２，…，ａｐ］Ｔを３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式を用いて前記演算部により計算し、前記ＡＲＭＡモデルにおけるＭＡ過程のｐ個のパラメータを行列要素とする行列ｂ＝［ｂ１，ｂ２，…，ｂｐ］Ｔを、前記計算により求めた行列ａ及び数１５を用いて、
【００３０】
【数１５】
2004302690.tif 000012
【００３１】
と演算部にて計算し、計算した行列ａ及び行列ｂに含まれる前記ＡＲＭＡモデルのパラメータ、並びに前記時系列の値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を演算部にて計算し、計算した新たな白色雑音を用いて新たな行列ａ及び行列ｂを演算部にて計算し、行列ａ及び行列ｂが収束するまで行列ａ及び行列ｂの計算を演算部にて繰り返すことを特徴とする。
【００３２】
本発明に係るパラメータ推定方法は、行列ａ＝［ａ１，ａ２，…，ａｐ］Ｔは、数１６に基づいて前記計算部により計算される
【００３３】
【数１６】
2004302690.tif 000013
【００３４】
ことを特徴とする。
【００３５】
本発明に係るデータ予測方法は、記憶部及び演算部を備えたコンピュータを用いて、複数のデータから構成される時系列から該時系列に加わる新たなデータを予測する方法において、時系列を記憶部に記憶し、前記パラメータ推定方法を用い、前記時系列をＡＲＭＡモデルに当てはめたときのパラメータを演算部にて計算し、計算したパラメータ、及び記憶している前記時系列の値を前記ＡＲＭＡモデルの式に代入して、前記時系列に加わる新たなデータを演算部にて計算することを特徴とする。
【００３６】
本発明に係るパラメータ推定装置は、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときのパラメータを推定する装置において、Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を記憶する手段と、平均値が０であるＮ個の値が時系列で並んだ白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成する手段と、生成した白色雑音の分散σｅ２を計算する手段と、前記ＡＲＭＡモデルにおけるＡＲ過程のｐ個のパラメータを行列要素とする行列ａ＝［ａ１，ａ２，…，ａｐ］Ｔを３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式を用いて前記演算部により計算する手段と、前記ＡＲＭＡモデルにおけるＭＡ過程のｐ個のパラメータを行列要素とする行列ｂ＝［ｂ１，ｂ２，…，ｂｐ］Ｔを、前記計算により求めた行列ａ及び数１７を用いて、
【００３７】
【数１７】
2004302690.tif 000014
【００３８】
と演算部にて計算する手段と、計算した行列ａ及び行列ｂに含まれる前記ＡＲＭＡモデルのパラメータ、並びに前記時系列の値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算する手段と、計算した新たな白色雑音を用いて新たな行列ａ及び行列ｂを計算する手段と、行列ａ及び行列ｂが収束するまで行列ａ及び行列ｂの計算を繰り返す手段とを備えることを特徴とする。
【００３９】
本発明に係るデータ予測装置は、複数のデータから構成される時系列から該時系列に加わる新たなデータを予測する装置において、時系列を記憶する手段と、前記パラメータ推定装置と、該パラメータ推定装置を用いて、前記時系列から、前記時系列をＡＲＭＡモデルに当てはめたときのパラメータを計算する手段と、計算したパラメータ、及び記憶している前記時系列の値を前記ＡＲＭＡモデルの式に代入して、前記時系列に加わる新たなデータを計算する手段とを備えることを特徴とする。
【００４０】
本発明に係るコンピュータプログラムは、コンピュータに、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときのパラメータを推定させるコンピュータプログラムにおいて、コンピュータに、平均値が０であるＮ個の値が時系列で並んだ白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成させる手順と、コンピュータに、生成した白色雑音の分散σｅ２を計算させる手順と、コンピュータに、記憶している、Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を用いて、前記ＡＲＭＡモデルにおけるＡＲ過程のｐ個のパラメータを行列要素とする行列ａ＝［ａ１，ａ２，…，ａｐ］Ｔを３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式により計算させる手順と、コンピュータに、記憶している、Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を用いて、前記ＡＲＭＡモデルにおけるＭＡ過程のｐ個のパラメータを行列要素とする行列ｂ＝［ｂ１，ｂ２，…，ｂｐ］Ｔを、前記計算により求めた行列ａ及び数１８により、
【００４１】
【数１８】
2004302690.tif 000015
【００４２】
と計算させる手順と、コンピュータに、計算した行列ａ及び行列ｂに含まれる前記ＡＲＭＡモデルのパラメータ、並びに前記時系列の値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算させる手順と、コンピュータに、計算した新たな白色雑音を用いて新たな行列ａ及び行列ｂを計算させる手順と、コンピュータに、行列ａ及び行列ｂが収束するまで行列ａ及び行列ｂの計算を繰り返させる手順とを含むことを特徴とする。
【００４３】
本発明に係るコンピュータプログラムは、コンピュータに、複数のデータから構成される時系列から該時系列に加わる新たなデータを予測させるコンピュータプログラムにおいて、コンピュータに、前記コンピュータプログラムを用いて、記憶している時系列から、該時系列をＡＲＭＡモデルに当てはめたときのパラメータを計算させる手順と、コンピュータに、計算したパラメータ、及び記憶している前記時系列の値を前記ＡＲＭＡモデルの式に代入して、前記時系列に加わる新たなデータを計算させる手順とを含むことを特徴とする。
【００４４】
本発明に係るコンピュータでの読取りが可能な記録媒体は、コンピュータに、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときのパラメータを推定させるコンピュータプログラムを記録してあるコンピュータでの読み取りが可能な記録媒体において、コンピュータに、平均値が０であるＮ個の値が時系列で並んだ白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成させる手順と、コンピュータに、生成した白色雑音の分散σｅ２を計算させる手順と、コンピュータに、記憶している、Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を用いて、前記ＡＲＭＡモデルにおけるＡＲ過程のｐ個のパラメータを行列要素とする行列ａ＝［ａ１，ａ２，…，ａｐ］Ｔを３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式により計算させる手順と、コンピュータに、記憶している、Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を用いて、前記ＡＲＭＡモデルにおけるＭＡ過程のｐ個のパラメータを行列要素とする行列ｂ＝［ｂ１，ｂ２，…，ｂｐ］Ｔを、前記計算により求めた行列ａ及び数１９により、
【００４５】
【数１９】
2004302690.tif 000016
【００４６】
と計算させる手順と、コンピュータに、計算した行列ａ及び行列ｂに含まれる前記ＡＲＭＡモデルのパラメータ、並びに前記時系列の値を前記ＡＲＭＡモデルの式に代入して、新たな白色雑音を計算させる手順と、コンピュータに、計算した新たな白色雑音を用いて新たな行列ａ及び行列ｂを計算させる手順と、コンピュータに、行列ａ及び行列ｂが収束するまで行列ａ及び行列ｂの計算を繰り返させる手順とを含むコンピュータプログラムを記録してあることを特徴とする。
【００４７】
本発明にあっては、３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式を用いることにより、ＡＲＭＡモデルのＭＡ項が正規性雑音の場合、ＡＲパラメータが不偏推定となる性質を利用する。つまり、（１）式のＡＲＭＡモデルのＭＡ項は（１６）式で表される。
【００４８】
【数２０】
2004302690.tif 000017
【００４９】
ここで、（１６）式のＭＡ項をＵｎとした場合、Ｕｎは有色ノイズである。しかしながら、ＵｎとＵｎ−１とは当然に相関を持つことから３次以上のキュムラントではＵｎをホワイトノイズとみなすことができる。そして、Ｕｎがホワイトノイズとみなせる場合は、３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式を用いて、まず第１に不偏推定値であるＡＲパラメータを推定することができる。
【００５０】
より詳細に説明すると、確率変数が正規分布に従う場合、３次以上のキュムラントは０となる。（１）式の両辺にＸｎ−１ ２を乗じて平均をとると｛ｅｎ｝が３次白色であるので（１７）式を得ることができる。
【００５１】
【数２１】
2004302690.tif 000018
【００５２】
ここで，ｃｌは（１８）式の如く３次のキュムラントである。
【００５３】
【数２２】
2004302690.tif 000019
【００５４】
（１７）式においてｌを１からｐまで動かすと、（１９）式の如く３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式を得ることができる。
【００５５】
【数２３】
2004302690.tif 000020
【００５６】
（１９）式は下記の（２０）式で表すことができ、これを用いて（２１）式を得る。
【００５７】
【数２４】
2004302690.tif 000021
【００５８】
（２１）式から、ＡＲパラメータ｛ａｉ  ｝を不偏推定値として求めることができる。続いて、ＭＡパラメータ｛ｂｊ  ｝について説明する。
【００５９】
本発明者は、ＡＲＭＡモデルの式の構成から導かれるａ  ｐｒｉｏｒｉ  ｋｎｏｗｌｅｄｇｅを導入することにより、従来の疑似線形回帰法の改良を試みた。（１）式より、ＡＲＭＡモデルでは、現在のデータの値は、過去のデータの値および過去の白色雑音の値に関連づけられている。しかし、現在のデータの値は、未来の白色雑音の値とは関連性がないため、データと該データよりも未来の白色雑音との相互相関関数は、理論的に、
Ｅ［ｘｎ  ｅｎ＋ｋ  ］＝０
となる。同様に、
Ｅ［ｅｎ  ｘｎ−ｋ  ］＝０
である。これにより、（５）式に含まれるＴ−ｋは理論的に０とすることができる。また、白色雑音｛ｅｎ  ｝の自己相関関数は、理論的に、
【００６０】
【数２５】
2004302690.tif 000022
【００６１】
となる。実際の場合は、Ｎの値が有限であるため、相互相関関数および自己相関関数を計算した値は、理論的な値とは違った値となる。理論的な値が得られる部分は、理論的な値と置き換えることにより、（５）式は、
【００６２】
【数２６】
2004302690.tif 000023
【００６３】
と表すことができる。単位行列Ｉを用いて、（２２）式は、更に、
【００６４】
【数２７】
2004302690.tif 000024
【００６５】
と表すことができる。（２３）式を、ベクトルａ及びｂについての連立方程式の形に書き直すと、
Ｒ・ａ−Ｔ・ｂ＝−ｒ        …（２４）
−ＴＴ  ・ａ＋σｅ２ｂ＝ｔ      …（２５）
となる。（２５）式より、
ｂ＝（ｔ＋ＴＴ  ・ａ）／σｅ２    …（２６）
が得られる。
【００６６】
本発明にあっては、演算部は（１９）式で表される３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式及び（２１）式を用いてまず、ＡＲパラメータであるベクトルａを計算する。さらに演算部は白色雑音｛ｅｎ  ｝を、乱数等を用いて生成し、生成した白色雑音と記憶部に記憶した時系列とから、（２６）式を用いてベクトルｂを計算する。演算部は、計算したベクトルａ及びｂから、（１）式を用いて新たな白色雑音｛ｅｎ  ｝を計算し、推定すべきパラメータであるベクトルｂが収束するまで計算を繰り返す。以上のようにベクトルａについては容易に求めることができ、またベクトルｂについては、Ｔが上三角行列であるため真値への収束が早く、計算時間を短縮することができるほか解のロバスト性を期待することが可能となる。
【００６７】
また、本発明にあっては従来より確実にパラメータを推定できる方法により推定したパラメータを用いてデータを予測することにより、従来ではデータの予測ができなかった時系列についてもデータを予測することができる。
【００６８】
【発明の実施の形態】
以下本発明を実施の形態を示す図面に基づいて詳述する。
図１は、本発明のデータ予測装置の構成を示すブロック図である。データ予測装置１は、本発明のパラメータ推定装置の機能を兼ね備えており、例えばパーソナルコンピュータを用いて構成されている。以下では、データ予測装置及びパラメータ推定装置をパーソナルコンピュータ１であるものとして説明する。
【００６９】
パーソナルコンピュータ１は、バス１７を介して演算部としてのＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）１１、キーボード、マウス等からなる入力部１３，液晶ディスプレイ等の表示部１４、ハードディスク等の記憶部１５，演算に伴って発生する一時的な情報を記憶するＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）１２，モデム等の通信部１６，ＣＤ−ＲＯＭドライブ等の外部記憶装置１８等が接続されている。ＣＰＵ１１は、バス１７を介して上述したハードウェア各部と接続されており、ＣＰＵ１１はそれらを制御すると共に、記憶部１５に格納された制御プログラム１５Ｐに従って、種々のソフトウェア的機能を実行する。制御プログラム１５Ｐは、以下に述べるソフトウェア処理をＣ言語等のプログラミング言語で記述されている。
【００７０】
また、ＣＤ−ＲＯＭ等の本発明に係る記録媒体２から本発明に係るコンピュータプログラム２０を外部記憶装置１８にて読み取り、読み取ったコンピュータプログラム２０を記憶部１５にインストールする形態であっても良い。この場合、記憶部１５にインストールされた制御プログラムは、ＣＰＵ１１の指示によりＲＡＭ１２にロードされ、ＣＰＵ１１は、ロードしたコンピュータプログラム２０に基づいて必要な処理を実行する。通信部１６には通信網Ｎ（以下、インターネットＮという）を介して、本発明に係るコンピュータプログラム２０を格納したＷＷＷ（Ｗｏｒｌｄ  Ｗｉｄｅ  Ｗｅｂ）サーバＳが接続されている。本発明にあっては、ＣＤ−ＲＯＭ等の記録媒体２によりインストールするほか、ＷＷＷサーバＳからコンピュータプログラム２０をダウンロードし、記憶部１５にインストールする形態であってもよい。
【００７１】
ＷＷＷサーバＳにはデータの予測の基礎となる時系列データが記憶されており、パーソナルコンピュータ１はＷＷＷサーバＳから配信される時系列データを利用して将来的なデータを予測する。以下では、取り扱うデータとして株価のデータに適用した場合について説明するが、これに限られるものではなく売上高の他、データは、入力部１３からオペレータにより入力されるデータまたは図示しないセンサから出力される機器の各種データ（温度データ、流量データ）等であっても良い。
【００７２】
図２は、パーソナルコンピュータ１が行うパラメータ推定の処理およびデータ予測の処理を示すフローチャートである。パーソナルコンピュータ１のＣＰＵ１１は、ＷＷＷサーバＳからデータをダウンロードする。ＣＰＵ１１はダウンロードしたＮ個のデータから構成されるデータの時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝をＲＡＭ１２に読み込む（ステップＳ２１）。ＣＰＵ１１は、コンピュータプログラム２０をＲＡＭ１２へロードする。まず、ＣＰＵ１１は（１９）式及び（２０）式に示されるベクトルｃ及びベクトルｃｓを計算する（ステップＳ２２）。具体的には、ＣＰＵ１１はＲＡＭ１２に読み込んだ時系列データ｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を（１８）式に代入することにより計算する。
【００７３】
続いて、ＣＰＵ１１は計算したベクトルｃ及びベクトルｃｓを（２１）式に代入することにより、ＡＲパラメータであるベクトルａを計算する（ステップＳ２３）。なお、（１８）式、（１９）式、（２０）式及び（２１）式等の計算に必要な式は記憶部１５に記憶されており、適宜ＲＡＭ１２に読み出されてＣＰＵ１１にて計算される。次に、ＣＰＵ１１は、Ｎ個の値から構成される平均値０の白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成する（ステップＳ２４）。このとき、ＣＰＵ１１は、乱数生成のルーチンを用い、時刻の値を初期値として乱数を発生させる等して、白色雑音を生成する。
【００７４】
ＣＰＵ１１は、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、生成した白色雑音｛ｅｎ  ｝の分散σｅ２を計算し（ステップＳ２５）、データの時系列｛ｘｎ  ｝及び白色雑音｛ｅｎ  ｝から、（６）式で定義された行列ｔと、（２３）式で定義された行列Ｔとを計算する（ステップＳ２６）。ＣＰＵ１１は、計算した行列ｔ，Ｔ及びステップＳ２３で計算したベクトルａ、ステップＳ２５で計算した白色雑音の分散σｅ２を、（２６）式に代入して、求めるべきＭＡパラメータであるベクトルｂを計算する（ステップＳ２７）。
【００７５】
ＣＰＵ１１は、次に、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、計算したベクトルａ及びｂが収束したか否かを判定する（ステップＳ２８）。ＣＰＵ１１は計算したベクトルａ及びｂが収束していないと判断した場合は（ステップＳ２８でＮＯ）、ベクトルａ及びｂとデータの時系列｛ｘｎ  ｝とを（１）式に代入して、新たな白色雑音｛ｅｎ  ｝を計算し（ステップＳ２９）、処理をステップＳ２５へ戻す。そして、ＣＰＵ１１は、ベクトルａ及びｂが収束するまでステップＳ２５〜Ｓ２９の処理を繰り返す。なお、（６）式、（２３）式、（３６）式及び（１）式等の計算に必要な式は記憶部１５に記憶されており、適宜ＲＡＭ１２に読み出されてＣＰＵ１１にて計算される。
【００７６】
ステップＳ２８においては、ＣＰＵ１１は、例えば、最新のベクトルａ，ｂと前回に計算したベクトルａ，ｂであるベクトルａ’，ｂ’とから、
Δ＝｜ａ−ａ’｜＋｜ｂ−ｂ’｜
を計算し、Δの値が所定値よりも小さい場合に、ベクトルａ及びｂが収束したと判定するようにすればよい。なお、所定値は予め記憶部１５に記憶されている。
【００７７】
ステップＳ２８において、ＣＰＵ１１はベクトルａ及びｂが収束していると判断した場合（ステップＳ２８でＹＥＳ）、ＣＰＵ１１は、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、計算したパラメータを用いて、ＡＲＭＡモデルにより時系列｛ｘｎ  ｝の予測データを計算し（ステップＳ２１０）、処理を終了する。以上の処理により、ＡＲＭＡモデルのパラメータが推定され、時系列の新たなデータが予測できる。
【００７８】
【発明の効果】
以上詳述した如く、本発明にあっては、３次キュムラントに対するＹｕｌｅ−Ｗａｌｋｅｒ方程式を用いてまず、ＡＲパラメータであるベクトルａを計算する。さらに白色雑音｛ｅｎ  ｝を、乱数等を用いて生成し、生成した白色雑音と記憶部に記憶した時系列とから、ベクトルｂを計算する。演算部は、計算したベクトルａ及びｂから、（１）式を用いて新たな白色雑音｛ｅｎ  ｝を計算し、推定すべきパラメータであるベクトルｂが収束するまで計算を繰り返す。従って、ベクトルａについては容易に求めることができ、またベクトルｂについては、Ｔが上三角行列であるため真値への収束が早く、計算時間を短縮することができるほか解のロバスト性を期待することが可能となる。
【００７９】
また、本発明にあっては従来より確実にパラメータを推定できる方法により推定したパラメータを用いてデータを予測することにより、従来ではデータの予測ができなかった時系列についてもデータを予測することができる等、本発明は優れた効果を奏し得る。
【図面の簡単な説明】
【図１】本発明のデータ予測装置の構成を示すブロック図である。
【図２】パーソナルコンピュータが行うパラメータ推定の処理およびデータ予測の処理を示すフローチャートである。
【符号の説明】
１  パーソナルコンピュータ（データ予測装置、パラメータ推定装置）
１１  ＣＰＵ（演算部）
１２  ＲＡＭ
１５  記憶部
１３  入力部
１４  表示部
１６  通信部
Ｓ  ＷＷＷサーバ
１５Ｐ  制御プログラム
２  記録媒体
２０  コンピュータプログラム
【図１】
2004302690.tif 000030
【図２】
2004302690.tif 000031