(57)【要約】【課題】従来よりも確実に時系列のパラメータを推定する方法、推定したパラメータよりデータを予測する方法、パラメータ推定装置、データ予測装置、コンピュータプログラム、及び記録媒体を提供する。【解決手段】データ予測装置１は、複数のデータからなる時系列を受け付け、白色雑音を生成し、ラグ０の白色雑音の自己相関関数を分散Ｓ０ ２  、ラグが０以外の白色雑音の自己相関関数を０、時系列と該時系列よりも未来の白色雑音の値との相互相関関数の値を０としてガウス−ニュートン法に必要な行列φ及びＹを計算し、時系列をＡＲＭＡモデルに当てはめたときのパラメータθの暫定値θ（ｈ）  から、行列ｆ（θ（ｈ）  ）及びεを計算し、ガウス−ニュートン法により次の暫定値θ（ｈ＋１）  を計算し、繰り返し計算によりθを決定し、θを用いて予測データを計算する。【選択図】        図１
2004287734.tif 000002
【特許請求の範囲】
【請求項１】
複数のデータからなる時系列、及び平均値が０で分散がＳ０２である白色雑音から、前記時系列の自己相関関数、前記白色雑音の自己相関関数、及び前記時系列と前記白色雑音との相互相関関数を用いて、前記時系列をＡＲＭＡモデルに当てはめたときのパラメータをガウス−ニュートン法により推定する方法において、
ラグが０の前記白色雑音の自己相関関数の値をＳ０ ２  とし、
ラグが０以外の前記白色雑音の自己相関関数の値を０とし、
前記時系列と前記白色雑音との相互相関関数のうち、前記時系列に含まれるデータと該データよりも未来の時点の白色雑音の値との相互相関関数の値を０として、
前記パラメータを計算することを特徴とするパラメータ推定方法。
【請求項２】
記憶部及び演算部を備えたコンピュータを用いて、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときの、ＡＲ過程のｐ個のパラメータ（ａ１  ，ａ２  ，…，ａｐ  ）及びＭＡ過程のｐ個のパラメータ（ｂ１  ，ｂ２  ，…，ｂｐ  ）を要素として含むパラメータθ＝［ａ１  ，ａ２  ，…，ａｐ  ，ｂ１  ，ｂ２  ，…，ｂｐ  ］Ｔ  を推定する方法において、
Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を記憶部に記憶し、
Ｒｋ  ＝Ｅ［ｘｎ  ｘｎ＋ｋ  ］で定義される｛ｘｎ  ｝の自己相関関数｛Ｒｋ  ：ｋ＝０，１，…，ｐ｝を演算部にて計算し、
平均値が０であるＮ個の値からなる白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を演算部にて生成し、
生成した白色雑音｛ｅｎ  ｝の分散Ｓ０ ２  を演算部にて計算し、
Ｔｋ  ＝Ｅ［ｅｎ  ｘｎ＋ｋ  ］で定義される｛ｅｎ  ｝と｛ｘｎ  ｝との相互相関関数｛Ｔｋ  ：ｋ＝０，１，…，ｐ｝を演算部にて計算し、
ｐ×ｐの単位行列Ｉを含む２ｐ×２ｐの行列φを、
【数１】
2004287734.tif 000028
と演算部にて生成し、
θの初期値を演算部にて生成してθの暫定値として記憶部に記憶し、
ｈ番目のθの暫定値θ（ｈ）  に対して、ε＝［ε１  ，ε２  ，…，εｐ  ，…，εｐ＋ｐ  ］Ｔ
を、
【数２】
2004287734.tif 000029
と演算部にて計算し、
Ｊ＝εＴ  εを演算部にて計算し、
計算したＪの値が所定の微少値より小さいか否かを演算部にて判定し、
Ｊの値が所定の微少値より小さいと判定された場合は、θの暫定値θ（ｈ）  がパラメータθの推定値であると演算部にて決定し、
Ｊの値が所定の微少値より小さいと判定されなかった場合は、（ｈ＋１）番目のθの暫定値θ（ｈ＋１）  を、
θ（ｈ＋１）  ＝θ（ｈ）  ＋φ−１ε
と演算部にて計算する
ことを特徴とするパラメータ推定方法。
【請求項３】
記憶部及び演算部を備えたコンピュータを用いて、複数のデータから構成される時系列から該時系列に加わる新たなデータを予測する方法において、
時系列を記憶部に記憶し、
請求項２に記載のパラメータ推定方法を用い、前記時系列をＡＲＭＡモデルに当てはめたときのパラメータを演算部にて計算し、
計算したパラメータ、及び記憶している前記時系列の値を前記ＡＲＭＡモデルの式に代入して、前記時系列に加わる新たなデータを演算部にて計算する
ことを特徴とするデータ予測方法。
【請求項４】
Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときの、ＡＲ過程のｐ個のパラメータ（ａ１  ，ａ２  ，…，ａｐ）及びＭＡ過程のｐ個のパラメータ（ｂ１  ，ｂ２  ，…，ｂｐ  ）を要素として含むパラメータθ＝［ａ１  ，ａ２  ，…，ａｐ  ，ｂ１  ，ｂ２  ，…，ｂｐ  ］Ｔ  を推定する装置において、
Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を記憶する手段と、
Ｒｋ  ＝Ｅ［ｘｎ  ｘｎ＋ｋ  ］で定義される｛ｘｎ  ｝の自己相関関数｛Ｒｋ  ：ｋ＝０，１，…，ｐ｝を計算する手段と、
平均値が０であるＮ個の値からなる白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成する手段と、
生成した白色雑音｛ｅｎ  ｝の分散Ｓ０２を計算する手段と、
Ｔｋ  ＝Ｅ［ｅｎ  ｘｎ＋ｋ  ］で定義される｛ｅｎ  ｝と｛ｘｎ  ｝との相互相関関数｛Ｔｋ  ：ｋ＝０，１，…，ｐ｝を計算する手段と、
ｐ×ｐの単位行列Ｉを含む２ｐ×２ｐの行列φを、
【数３】
2004287734.tif 000030
と生成する手段と、
θの初期値を生成してθの最初の暫定値とする手段と、
ｈ番目のθの暫定値θ（ｈ）  に対して、ε＝［ε１  ，ε２  ，…，εｐ  ，…，εｐ＋ｐ  ］Ｔ
を、
【数４】
2004287734.tif 000031
と計算する手段と、
Ｊ＝εＴ  εを計算する手段と、
計算したＪの値が所定の微少値より小さいか否かを判定する手段と、
Ｊの値が所定の微少値より小さいと判定された場合は、θの暫定値θ（ｈ）  がパラメータθの推定値であると決定する手段と、
Ｊの値が所定の微少値より小さいと判定されなかった場合は、（ｈ＋１）番目のθの暫定値θ（ｈ＋１）  を、
θ（ｈ＋１）  ＝θ（ｈ）  ＋φ−１ε
と計算する手段と
を備えることを特徴とするパラメータ推定装置。
【請求項５】
複数のデータから構成される時系列から該時系列に加わる新たなデータを予測する装置において、
時系列を記憶する手段と、
請求項４に記載のパラメータ推定装置と、
該パラメータ推定装置を用いて、前記時系列から、前記時系列をＡＲＭＡモデルに当てはめたときのパラメータを計算する手段と、
計算したパラメータ、及び記憶している前記時系列の値を前記ＡＲＭＡモデルの式に代入して、前記時系列に加わる新たなデータを計算する手段と
を備えることを特徴とするデータ予測装置。
【請求項６】
コンピュータに、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときの、ＡＲ過程のｐ個のパラメータ（ａ１  ，ａ２  ，…，ａｐ  ）及びＭＡ過程のｐ個のパラメータ（ｂ１  ，ｂ２  ，…，ｂｐ  ）を要素として含むパラメータθ＝［ａ１  ，ａ２  ，…，ａｐ  ，ｂ１  ，ｂ２  ，…，ｂｐ  ］Ｔ  を推定させるコンピュータプログラムにおいて、
コンピュータに、Ｎ個のデータから構成された記憶している時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝から、Ｒｋ  ＝Ｅ［ｘｎ  ｘｎ＋ｋ  ］で定義される｛ｘｎ  ｝の自己相関関数｛Ｒｋ  ：ｋ＝０，１，…，ｐ｝を計算させる手順と、
コンピュータに、平均値が０であるＮ個の値からなる白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成させる手順と、
コンピュータに、生成した白色雑音｛ｅｎ  ｝の分散Ｓ０２を計算させる手順と、
コンピュータに、Ｔｋ  ＝Ｅ［ｅｎ  ｘｎ＋ｋ  ］で定義される｛ｅｎ  ｝と｛ｘｎ  ｝との相互相関関数｛Ｔｋ  ：ｋ＝０，１，…，ｐ｝を計算させる手順と、
コンピュータに、ｐ×ｐの単位行列Ｉを含む２ｐ×２ｐの行列φを、
【数５】
2004287734.tif 000032
と生成させる手順と、
コンピュータに、θの初期値を生成してθの最初の暫定値とさせる手順と、
コンピュータに、ｈ番目のθの暫定値θ（ｈ）  に対して、ε＝［ε１  ，ε２  ，…，εｐ  ，…，εｐ＋ｐ  ］Ｔ  を、
【数６】
2004287734.tif 000033
と計算させる手順と、
コンピュータに、Ｊ＝εＴ  εを計算させる手順と、
コンピュータに、計算したＪの値が所定の微少値より小さいか否かを判定させる手順と、
Ｊの値が所定の微少値より小さいと判定された場合は、コンピュータに、θの暫定値θ（ｈ）  がパラメータθの推定値であると決定させる手順と、
Ｊの値が所定の微少値より小さいと判定されなかった場合は、コンピュータに、（ｈ＋１）番目のθの暫定値θ（ｈ＋１）  を、
θ（ｈ＋１）  ＝θ（ｈ）  ＋φ−１ε
と計算させる手順と
を含むことを特徴とするコンピュータプログラム。
【請求項７】
コンピュータに、複数のデータから構成される時系列から該時系列に加わる新たなデータを予測させるコンピュータプログラムにおいて、
コンピュータに、請求項６に記載のコンピュータプログラムを用いて、記憶している時系列から、該時系列をＡＲＭＡモデルに当てはめたときのパラメータを計算させる手順と、
コンピュータに、計算したパラメータ、及び記憶している前記時系列の値を前記ＡＲＭＡモデルの式に代入して、前記時系列に加わる新たなデータを計算させる手順と
を含むことを特徴とするコンピュータプログラム。
【請求項８】
コンピュータに、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときの、ＡＲ過程のｐ個のパラメータ（ａ１  ，ａ２  ，…，ａｐ  ）及びＭＡ過程のｐ個のパラメータ（ｂ１  ，ｂ２  ，…，ｂｐ  ）を要素として含むパラメータθ＝［ａ１  ，ａ２  ，…，ａｐ  ，  ｂ１  ，ｂ２  ，…，ｂｐ  ］Ｔ  を推定させるコンピュータプログラムを記録してあるコンピュータでの読み取りが可能な記録媒体において、
コンピュータに、Ｎ個のデータから構成された記憶している時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝から、Ｒｋ  ＝Ｅ［ｘｎ  ｘｎ＋ｋ  ］で定義される｛ｘｎ  ｝の自己相関関数｛Ｒｋ  ：ｋ＝０，１，…，ｐ｝を計算させる手順と、
コンピュータに、平均値が０であるＮ個の値からなる白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成させる手順と、
コンピュータに、生成した白色雑音｛ｅｎ  ｝の分散Ｓ０２を計算させる手順と、
コンピュータに、Ｔｋ  ＝Ｅ［ｅｎ  ｘｎ＋ｋ  ］で定義される｛ｅｎ  ｝と｛ｘｎ  ｝との相互相関関数｛Ｔｋ  ：ｋ＝０，１，…，ｐ｝を計算させる手順と、
コンピュータに、ｐ×ｐの単位行列Ｉを含む２ｐ×２ｐの行列φを、
【数７】
2004287734.tif 000034
と生成させる手順と、
コンピュータに、θの初期値を生成してθの最初の暫定値とさせる手順と、
コンピュータに、ｈ番目のθの暫定値θ（ｈ）  に対して、ε＝［ε１  ，ε２  ，…，εｐ  ，…，εｐ＋ｐ  ］Ｔ  を、
【数８】
2004287734.tif 000035
と計算させる手順と、
コンピュータに、Ｊ＝εＴ  εを計算させる手順と、
コンピュータに、計算したＪの値が所定の微少値より小さいか否かを判定させる手順と、
Ｊの値が所定の微少値より小さいと判定された場合は、コンピュータに、θの暫定値θ（ｈ）  がパラメータθの推定値であると決定させる手順と、
Ｊの値が所定の微少値より小さいと判定されなかった場合は、コンピュータに、（ｈ＋１）番目のθの暫定値θ（ｈ＋１）  を、
θ（ｈ＋１）  ＝θ（ｈ）  ＋φ−１ε
と計算させる手順と
を含むコンピュータプログラムを記録してあることを特徴とするコンピュータでの読み取りが可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、データの時系列からパラメータを推定する方法、推定されたパラメータを用いて新たなデータを予測する方法、パラメータ推定装置、データ予測装置、コンピュータをパラメータ推定装置またはデータ予測装置として実現するためのコンピュータプログラム、及び記録媒体に関する。
【０００２】
【従来の技術】
蓄積されたデータの時系列を解析し、新たなデータの値を予想する方法は、従来より様々な方法が提案されている。株価または販売量などのデータから構成される一般の経済的な時系列は、白色ノイズによって励起された系の出力とみることができる。データの時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ，…｝を、定常エルゴード的正規過程ｘ（ｔ）の標本時系列であるとしたときには、（ｐ，ｑ）次のＡＲＭＡモデルでは、｛ｘｎ  ｝の要素間の関係は、以下の式で表される。
【０００３】
【数９】
2004287734.tif 000003
【０００４】
ここで、｛ｅｎ  ｝は、平均値０及び分散σｅ ２  の正規性雑音であり、｛ａｉ  ：ｉ＝１，２，…，ｐ｝及び｛ｂｊ  ：ｊ＝１，２，…，ｑ｝は、
Ａ（Ｚ−１）＝１＋ａ１  Ｚ−１＋…＋ａｐ  Ｚ−ｐ
Ｂ（Ｚ−１）＝１＋ｂ１  Ｚ−１＋…＋ｂｑ  Ｚ−ｑ
とおいたときに、Ａ（Ｚ−１）及びＢ（Ｚ−１）が規約であり、更に定常条件、可逆条件および強正実条件を満たすものとする。
【０００５】
標本の時系列｛ｘｎ  ｝から（１）式のパラメータ｛ａｉ  ｝及び｛ｂｊ  ｝を推定するための方法として、非特許文献１に開示されたガウス−ニュートン法がある。説明変数ｃ１  ，…，ｃｍ  と目的変数ｙとに対して、測定値の組（ｃ１ｉ，…，ｃｍｉ；ｙｉ  ），ｉ＝１，…，Ｌが与えられているとすると、
2004287734.tif 000004
とする。ここで、εｉ  は、平均値０及び分散σｉ ２  のガウス分布に従い、かつ互いに独立であるとする。更に、
ｗｉ  ＝１／σｉ ２          ｉ＝１，…，Ｌ
の重みｗｉ  を考え、評価関数Ｓ（θ）を、
【０００６】
【数１０】
2004287734.tif 000005
【０００７】
と表す。Ｓ（θ）を最小にするθを求めることにより、θの最尤推定値を求めることができる。θを（θ１  ，…，θｍ  ）として、Ｓ（θ）のグラジエント∇Ｓ（θ）及びヘッシアン∇２  Ｓ（θ）を、
【０００８】
【数１１】
2004287734.tif 000006
【０００９】
と定義する。ガウス−ニュートン法によるアルゴリズムでは、繰り返し計算によってθを計算し、ｈ＋１番目のθであるθ（ｈ＋１）  は、
θ（ｈ＋１）  ＝θ（ｈ）  −［∇２Ｓ（θ（ｈ）  ）］−１∇Ｓ（θ（ｈ）  ）  …（５）
で得られる。ここで、
【００１０】
【数１２】
2004287734.tif 000007
【００１１】
とおくと、（５）式は、（３）式および（４）式より、
（φＴ  Ｗφ）Δθ（ｈ）  ＝φＴ  Ｗ［Ｙ−ｆ（θ（ｈ）  ）］  …（７）
と表される。
【００１２】
次に、ガウス−ニュートン法を、（１）式のパラメータを求める問題に適用する。簡単にするためにｐ＝ｑとし、ベクトルθを次の様に定義する。
θ＝［ａ１  ，ａ２  ，…，ａｐ  ，ｂ１  ，ｂ２  ，…，ｂｐ  ］Ｔ      …（８）
ラグがｋである｛ｘｎ  ｝の自己相関関数をＲｋ  とすると、
Ｒｋ  ＝Ｅ［ｘｎ  ｘｎ＋ｋ  ］
である。ここで、Ｅ［ｄｎ  ］は、データ列｛ｄｎ  ：ｎ＝１，２，…｝の期待値である。また、
Ｒ−ｋ＝Ｒｋ
である。ラグがｋである｛ｅｎ  ｝の自己相関関数をＳｋとすると、同様に、
Ｓｋ  ＝Ｅ［ｅｎ  ｅｎ＋ｋ  ］，Ｓ−ｋ＝Ｓｋ
である。また、｛ｘｎ  ｝と｛ｅｎ  ｝との相互相関関数は、
Ｔｋ  ＝Ｅ［ｅｎ  ｘｎ＋ｋ  ］
Ｔ−ｋ＝Ｅ［ｅｎ  ｘｎ−ｋ  ］
とおくと、
Ｅ［ｘｎ  ｅｎ＋ｋ  ］＝Ｅ［ｅｎ  ｘｎ−ｋ  ］＝Ｔ−ｋ
Ｅ［ｘｎ  ｅｎ−ｋ  ］＝Ｅ［ｅｎ  ｘｎ＋ｋ  ］＝Ｔｋ
となる。（１）式の両辺にｅｎ−１  ，…，ｅｎ−ｐ  を掛けてＥ［］の操作を行うと、
Ｔ０  ＋ａ１  Ｔ−１  ＋…＋ａｐ  Ｔ−ｐ＝Ｓ０  ＋ｂ１  Ｓ１  ＋…＋ｂｐ  Ｓｐ      …（９−０）
Ｔ１  ＋ａ１  Ｔ０  ＋…＋ａｐ  Ｔ−（ｐ−１）  ＝Ｓ１  ＋ｂ１  Ｓ０  ＋…＋ｂｐ  Ｓｐ−１      …（９−１）
Ｔ２  ＋ａ１  Ｔ１  ＋…＋ａｐ  Ｔ−（ｐ−２）  ＝Ｓ２  ＋ｂ１  Ｓ１  ＋…＋ｂｐ  Ｓｐ−２      …（９−２）
…
Ｔｐ  ＋ａ１  Ｔｐ−１  ＋…＋ａｐ  Ｔ０  ＝Ｓｐ  ＋ｂ１  Ｓｐ−１  ＋…＋ｂｐ  Ｓ０      …（９−ｐ）
が得られる。同様に、ｘｎ  ，ｘｎ−１  ，…，ｘｎ−ｐ  を（１）式の両辺に掛けてＥ［］の操作を行うと、
Ｒ０  ＋ａ１  Ｒ１  ＋…＋ａｐ  Ｒｐ  ＝Ｔ０  ＋ｂ１  Ｔ１  ＋…＋ｂｐ  Ｔｐ      …（１０−０）
Ｒ１  ＋ａ１  Ｒ０  ＋…＋ａｐ  Ｒｐ−１  ＝Ｔ−１＋ｂ１  Ｔ０  ＋…＋ｂｐ  Ｔｐ−１      …（１０−１）
Ｒ２  ＋ａ１  Ｒ１  ＋…＋ａｐ  Ｒｐ−２  ＝Ｔ−２＋ｂ１  Ｔ−１＋…＋ｂｐ  Ｔｐ−２      …（１０−２）
…
Ｒｐ  ＋ａ１  Ｒｐ−１  ＋…＋ａｐ  Ｒ０  ＝Ｔ−ｐ＋ｂ１  Ｔ−（ｐ−１）  ＋…＋ｂｐ  Ｔ０      …（１０−ｐ）
が得られる。（９−１）〜（９−ｐ）（１０−１）〜（１０−ｐ）式に含まれるパラメータの値を推定値とし、各パラメータの推定値から構成される各式の右辺および左辺の間の残差をεｉ  とすると、
Ｔ１  ＋ａ１  Ｔ０  ＋…＋ａｐ  Ｔ−（ｐ−１）  −（Ｓ１  ＋ｂ１  Ｓ０  ＋…＋ｂｐ  Ｓｐ−１  ）＝ε１
Ｔ２  ＋ａ１  Ｔ１  ＋…＋ａｐ  Ｔ−（ｐ−２）  −（Ｓ２  ＋ｂ１  Ｓ１  ＋…＋ｂｐ  Ｓｐ−２  ）＝ε２
…
Ｔｐ  ＋ａ１  Ｔｐ−１  ＋…＋ａｐ  Ｔ０  −（Ｓｐ  ＋ｂ１  Ｓｐ−１  ＋…＋ｂｐ  Ｓ０  ）＝εｐ
−（Ｒ１  ＋ａ１  Ｒ０  ＋…＋ａｐ  Ｒｐ−１  ）＋（Ｔ−１＋ｂ１  Ｔ０  ＋…＋ｂｐ  Ｔｐ−１  ）＝εｐ＋１
−（Ｒ２  ＋ａ１  Ｒ１  ＋…＋ａｐ  Ｒｐ−２  ）＋（Ｔ−２＋ｂ１  Ｔ−１＋…＋ｂｐ  Ｔｐ−２  ）＝εｐ＋２
…
−（Ｒｐ  ＋ａ１  Ｒｐ−１  ＋…＋ａｐ  Ｒ０  ）＋（Ｔ−ｐ＋ｂ１  Ｔ−（ｐ−１）  ＋…＋ｂｐ  Ｔ０  ）＝εｐ＋ｐ
が得られる。この式を更に変形すると、
【００１３】
【数１３】
2004287734.tif 000008
【００１４】
が得られる。（２）式の形式に合わせるため、
Ｙ＝［Ｔ１  ，Ｔ２  ，…，Ｔｐ  ，−Ｒ１  ，−Ｒ２  ，…，−Ｒｐ  ］Ｔ      …（１２）
ε＝［ε１  ，ε２  ，…，εｐ  ，εｐ＋１  ，εｐ＋２  ，…，εｐ＋ｐ  ］Ｔ    …（１３）
とし、更に（８）式のθに対して、
【００１５】
【数１４】
2004287734.tif 000009
【００１６】
とすると、（１１）式は、
Ｙ＝ｆ（θ）＋ε  …（１５）
となる。（１５）式のεを最小にするように、ガウス−ニュートン法を適用してθの推定値を求めれば良い。（１４）式を（６）式に適用して、このときのφは、
【００１７】
【数１５】
2004287734.tif 000010
【００１８】
となる。ここで、ｐ×ｐの行列Ｓ、τ及びＲを、
【００１９】
【数１６】
2004287734.tif 000011
【００２０】
と定義すると、φは、
【００２１】
【数１７】
2004287734.tif 000012
【００２２】
となる。φは２ｐ×２ｐの正方行列であるため、（７）式は、両辺のφＴ  Ｗが消えて、
θ（ｈ＋１）  ＝θ（ｈ）  ＋φ−１ε（ｈ）      …（１８）
となる。また、（１５）式より、
ε（ｈ）  ＝Ｙ−ｆ（θ（ｈ）  ）  …（１９）
である。
【００２３】
観測データ｛ｘｎ  ｝が得られた場合、｛ｘｎ  ｝の自己相関関数｛Ｒｋ  ｝を計算し、白色雑音｛ｅｎ  ｝を生成し、｛ｅｎ  ｝の自己相関関数｛Ｓｋ  ｝、及び｛ｘｎ  ｝と｛ｅｎ  ｝との相互相関関数｛Ｔｋ  ｝を計算し、パラメータの初期値θ（０）  を生成し、（１９）式を用いてε（ｈ）  を計算し、（１８）式を用いてθ（ｈ＋１）  を計算し、εが所定の微少値よりも小さくなるまで計算を繰り返して、パラメータθの推定値を得ることができる。
【００２４】
【非特許文献１】
得丸英勝、添田喬、中溝高好、秋月影雄「係数・測定−ランダムデータ処理の理論と応用」培風館  １９８２
【００２５】
【発明が解決しようとする課題】
従来の方法では、問題の性質によっては、パラメータθが収束せず、パラメータの推定ができない場合があるという問題がある。また、２ｐ行２ｐ列の行列計算を行うため、計算時間が長いという問題があった。
【００２６】
本発明は、以上の問題を解決するためになされたものであって、その目的とするところは、計算に用いる数値のうち理論的に０となるはずの数値を０とおいて計算を行うことにより、パラメータの値の収束を容易にし、同時に計算時間の短縮を図ることができるパラメータ推定方法、パラメータ推定装置、コンピュータを該パラメータ推定装置として実現するためのコンピュータプログラム、及び記録媒体を提供することにある。
【００２７】
また、本発明の他の目的とするところは、前記パラメータ推定方法により推定したパラメータを用いてデータの推定を行うことにより、従来の手法ではデータの予測ができなかった時系列についてもデータの予測が可能となるデータ予測方法、データ予測装置、コンピュータを該データ予測装置として実現するためのコンピュータプログラム、及び記録媒体を提供することにある。
【００２８】
【課題を解決するための手段】
第１発明に係るパラメータ推定方法は、複数のデータからなる時系列、及び平均値が０で分散がＳ０２である白色雑音から、前記時系列の自己相関関数、前記白色雑音の自己相関関数、及び前記時系列と前記白色雑音との相互相関関数を用いて、前記時系列をＡＲＭＡモデルに当てはめたときのパラメータをガウス−ニュートン法により推定する方法において、ラグが０の前記白色雑音の自己相関関数の値をＳ０２とし、ラグが０以外の前記白色雑音の自己相関関数の値を０とし、前記時系列と前記白色雑音との相互相関関数のうち、前記時系列に含まれるデータと該データよりも未来の時点の白色雑音の値との相互相関関数の値を０として、前記パラメータを計算することを特徴とする。
【００２９】
第２発明に係るパラメータ推定方法は、記憶部及び演算部を備えたコンピュータを用いて、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときの、ＡＲ過程のｐ個のパラメータ（ａ１  ，ａ２  ，…，ａｐ  ）及びＭＡ過程のｐ個のパラメータ（ｂ１  ，ｂ２  ，…，ｂｐ  ）を要素として含むパラメータθ＝［ａ１  ，ａ２  ，…，ａｐ  ，ｂ１  ，ｂ２  ，…，ｂｐ  ］Ｔ  を推定する方法において、Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を記憶部に記憶し、Ｒｋ  ＝Ｅ［ｘｎ  ｘｎ＋ｋ  ］で定義される｛ｘｎ  ｝の自己相関関数｛Ｒｋ  ：ｋ＝０，１，…，ｐ｝を演算部にて計算し、平均値が０であるＮ個の値からなる白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を演算部にて生成し、生成した白色雑音｛ｅｎ  ｝の分散Ｓ０２を演算部にて計算し、Ｔｋ  ＝Ｅ［ｅｎ  ｘｎ＋ｋ  ］で定義される｛ｅｎ  ｝と｛ｘｎ  ｝との相互相関関数｛Ｔｋ  ：ｋ＝０，１，…，ｐ｝を演算部にて計算し、ｐ×ｐの単位行列Ｉを含む２ｐ×２ｐの行列φを、
【００３０】
【数１８】
2004287734.tif 000013
【００３１】
と演算部にて生成し、θの初期値を演算部にて生成してθの暫定値として記憶部に記憶し、ｈ番目のθの暫定値θ（ｈ）  に対して、ε＝［ε１  ，ε２  ，…，εｐ  ，…，εｐ＋ｐ］Ｔ  を、
【００３２】
【数１９】
2004287734.tif 000014
【００３３】
と演算部にて計算し、Ｊ＝εＴεを演算部にて計算し、計算したＪの値が所定の微少値より小さいか否かを演算部にて判定し、Ｊの値が所定の微少値より小さいと判定された場合は、θの暫定値θ（ｈ）  がパラメータθの推定値であると演算部にて決定し、Ｊの値が所定の微少値より小さいと判定されなかった場合は、（ｈ＋１）番目のθの暫定値θ（ｈ＋１）  を、
θ（ｈ＋１）  ＝θ（ｈ）  ＋φ−１ε
と演算部にて計算することを特徴とする。
【００３４】
第３発明に係るデータ予測方法は、記憶部及び演算部を備えたコンピュータを用いて、複数のデータから構成される時系列から該時系列に加わる新たなデータを予測する方法において、時系列を記憶部に記憶し、第２発明に係るパラメータ推定方法を用い、前記時系列をＡＲＭＡモデルに当てはめたときのパラメータを演算部にて計算し、計算したパラメータ、及び記憶している前記時系列の値を前記ＡＲＭＡモデルの式に代入して、前記時系列に加わる新たなデータを演算部にて計算することを特徴とする。
【００３５】
第４発明に係るパラメータ推定装置は、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときの、ＡＲ過程のｐ個のパラメータ（ａ１  ，ａ２  ，…，ａｐ  ）及びＭＡ過程のｐ個のパラメータ（ｂ１  ，ｂ２  ，…，ｂｐ  ）を要素として含むパラメータθ＝［ａ１  ，ａ２  ，…，ａｐ  ，ｂ１  ，ｂ２  ，…，ｂｐ  ］Ｔを推定する装置において、Ｎ個のデータから構成される時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝を記憶する手段と、Ｒｋ  ＝Ｅ［ｘｎ  ｘｎ＋ｋ  ］で定義される｛ｘｎ  ｝の自己相関関数｛Ｒｋ  ：ｋ＝０，１，…，ｐ｝を計算する手段と、平均値が０であるＮ個の値からなる白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成する手段と、生成した白色雑音｛ｅｎ  ｝の分散Ｓ０２を計算する手段と、Ｔｋ  ＝Ｅ［ｅｎ  ｘｎ＋ｋ  ］で定義される｛ｅｎ  ｝と｛ｘｎ  ｝との相互相関関数｛Ｔｋ  ：ｋ＝０，１，…，ｐ｝を計算する手段と、ｐ×ｐの単位行列Ｉを含む２ｐ×２ｐの行列φを、
【００３６】
【数２０】
2004287734.tif 000015
【００３７】
と生成する手段と、θの初期値を生成してθの最初の暫定値とする手段と、ｈ番目のθの暫定値θ（ｈ）  に対して、ε＝［ε１  ，ε２  ，…，εｐ  ，…，εｐ＋ｐ  ］Ｔ  を、
【００３８】
【数２１】
2004287734.tif 000016
【００３９】
と計算する手段と、Ｊ＝εＴ  εを計算する手段と、計算したＪの値が所定の微少値より小さいか否かを判定する手段と、Ｊの値が所定の微少値より小さいと判定された場合は、θの暫定値θ（ｈ）  がパラメータθの推定値であると決定する手段と、Ｊの値が所定の微少値より小さいと判定されなかった場合は、（ｈ＋１）番目のθの暫定値θ（ｈ＋１）  を、
θ（ｈ＋１）  ＝θ（ｈ）  ＋φ−１ε
と計算する手段とを備えることを特徴とする。
【００４０】
第５発明に係るデータ予測装置は、複数のデータから構成される時系列から該時系列に加わる新たなデータを予測する装置において、時系列を記憶する手段と、第４発明に係るパラメータ推定装置と、該パラメータ推定装置を用いて、前記時系列から、前記時系列をＡＲＭＡモデルに当てはめたときのパラメータを計算する手段と、計算したパラメータ、及び記憶している前記時系列の値を前記ＡＲＭＡモデルの式に代入して、前記時系列に加わる新たなデータを計算する手段とを備えることを特徴とする。
【００４１】
第６発明に係るコンピュータプログラムは、コンピュータに、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときの、ＡＲ過程のｐ個のパラメータ（ａ１  ，ａ２  ，…，ａｐ  ）及びＭＡ過程のｐ個のパラメータ（ｂ１  ，ｂ２，…，ｂｐ  ）を要素として含むパラメータθ＝［ａ１  ，ａ２  ，…，ａｐ  ，ｂ１  ，ｂ２  ，…，ｂｐ  ］Ｔ  を推定させるコンピュータプログラムにおいて、コンピュータに、Ｎ個のデータから構成された記憶している時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝から、Ｒｋ  ＝Ｅ［ｘｎ  ｘｎ＋ｋ  ］で定義される｛ｘｎ  ｝の自己相関関数｛Ｒｋ  ：ｋ＝０，１，…，ｐ｝を計算させる手順と、コンピュータに、平均値が０であるＮ個の値からなる白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成させる手順と、コンピュータに、生成した白色雑音｛ｅｎ  ｝の分散Ｓ０ ２  を計算させる手順と、コンピュータに、Ｔｋ  ＝Ｅ［ｅｎ  ｘｎ＋ｋ  ］で定義される｛ｅｎ  ｝と｛ｘｎ  ｝との相互相関関数｛Ｔｋ  ：ｋ＝０，１，…，ｐ｝を計算させる手順と、コンピュータに、ｐ×ｐの単位行列Ｉを含む２ｐ×２ｐの行列φを、
【００４２】
【数２２】
2004287734.tif 000017
【００４３】
と生成させる手順と、コンピュータに、θの初期値を生成してθの最初の暫定値とさせる手順と、コンピュータに、ｈ番目のθの暫定値θ（ｈ）  に対して、ε＝［ε１  ，ε２  ，…，εｐ  ，…，εｐ＋ｐ  ］Ｔ  を、
【００４４】
【数２３】
2004287734.tif 000018
【００４５】
と計算させる手順と、コンピュータに、Ｊ＝εＴ  εを計算させる手順と、コンピュータに、計算したＪの値が所定の微少値より小さいか否かを判定させる手順と、Ｊの値が所定の微少値より小さいと判定された場合は、コンピュータに、θの暫定値θ（ｈ）  がパラメータθの推定値であると決定させる手順と、Ｊの値が所定の微少値より小さいと判定されなかった場合は、コンピュータに、（ｈ＋１）番目のθの暫定値θ（ｈ＋１）  を、
θ（ｈ＋１）  ＝θ（ｈ）  ＋φ−１ε
と計算させる手順とを含むことを特徴とする。
【００４６】
第７発明に係るコンピュータプログラムは、コンピュータに、複数のデータから構成される時系列から該時系列に加わる新たなデータを予測させるコンピュータプログラムにおいて、コンピュータに、第６発明に係るコンピュータプログラムを用いて、記憶している時系列から、該時系列をＡＲＭＡモデルに当てはめたときのパラメータを計算させる手順と、コンピュータに、計算したパラメータ、及び記憶している前記時系列の値を前記ＡＲＭＡモデルの式に代入して、前記時系列に加わる新たなデータを計算させる手順とを含むことを特徴とする。
【００４７】
第８発明に係るコンピュータでの読み取りが可能な記録媒体は、コンピュータに、Ｎ個のデータから構成される時系列を（ｐ，ｐ）次のＡＲＭＡモデルに当てはめたときの、ＡＲ過程のｐ個のパラメータ（ａ１  ，ａ２  ，…，ａｐ  ）及びＭＡ過程のｐ個のパラメータ（ｂ１  ，ｂ２  ，…，ｂｐ）を要素として含むパラメータθ＝［ａ１  ，ａ２  ，…，ａｐ  ，ｂ１  ，ｂ２  ，…，ｂｐ  ］Ｔ  を推定させるコンピュータプログラムを記録してあるコンピュータでの読み取りが可能な記録媒体において、コンピュータに、Ｎ個のデータから構成された記憶している時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝から、Ｒｋ  ＝Ｅ［ｘｎ  ｘｎ＋ｋ  ］で定義される｛ｘｎ  ｝の自己相関関数｛Ｒｋ  ：ｋ＝０，１，…，ｐ｝を計算させる手順と、コンピュータに、平均値が０であるＮ個の値からなる白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成させる手順と、コンピュータに、生成した白色雑音｛ｅｎ  ｝の分散Ｓ０ ２  を計算させる手順と、コンピュータに、Ｔｋ  ＝Ｅ［ｅｎ  ｘｎ＋ｋ  ］で定義される｛ｅｎ  ｝と｛ｘｎ  ｝との相互相関関数｛Ｔｋ  ：ｋ＝０，１，…，ｐ｝を計算させる手順と、コンピュータに、ｐ×ｐの単位行列Ｉを含む２ｐ×２ｐの行列φを、
【００４８】
【数２４】
2004287734.tif 000019
【００４９】
と生成させる手順と、コンピュータに、θの初期値を生成してθの最初の暫定値とさせる手順と、コンピュータに、ｈ番目のθの暫定値θ（ｈ）  に対して、ε＝［ε１  ，ε２  ，…，εｐ  ，…，εｐ＋ｐ  ］Ｔ  を、
【００５０】
【数２５】
2004287734.tif 000020
【００５１】
と計算させる手順と、コンピュータに、Ｊ＝εＴ  εを計算させる手順と、コンピュータに、計算したＪの値が所定の微少値より小さいか否かを判定させる手順と、Ｊの値が所定の微少値より小さいと判定された場合は、コンピュータに、θの暫定値θ（ｈ）  がパラメータθの推定値であると決定させる手順と、Ｊの値が所定の微少値より小さいと判定されなかった場合は、コンピュータに、（ｈ＋１）番目のθの暫定値θ（ｈ＋１）  を、
θ（ｈ＋１）  ＝θ（ｈ）  ＋φ−１ε
と計算させる手順とを含むコンピュータプログラムを記録してあることを特徴とする。
【００５２】
本発明者は、ＡＲＭＡモデルの式の構成から導かれるａ  ｐｒｉｏｒｉ  ｋｎｏｗｌｅｄｇｅを導入することにより、従来の疑似線形回帰法の改良を試みた。（１）式より、ＡＲＭＡモデルでは、現在のデータの値は、過去のデータの値および過去の白色雑音の値に関連づけられている。しかし、現在のデータの値は、未来の白色雑音の値とは関連性がないため、データと該データよりも未来の白色雑音との相互相関関数は、理論的に、
Ｅ［ｘｎ  ｅｎ＋ｋ  ］＝０
となる。同様に、
Ｅ［ｅｎ  ｘｎ−ｋ  ］＝０
である。これにより、（１１）式に含まれるＴ−ｋは理論的に０とすることができる。また、白色雑音｛ｅｎ  ｝の自己相関関数は、理論的に、
【００５３】
【数２６】
2004287734.tif 000021
【００５４】
となる。実際の場合は、Ｎの値が有限であるため、相互相関関数および自己相関関数を計算した値は、理論的な値とは違った値となる。理論的な値が得られる部分は、理論的な値と置き換え、σｅ ２  ＝Ｓ０  とすると、（９−０）〜（９−ｐ）式は、
Ｔ０                      ＝Ｓ０      …（２０−０）
Ｔ１  ＋ａ１  Ｔ０                        ＝ｂ１  Ｓ０      …（２０−１）
Ｔ２  ＋ａ１  Ｔ１  ＋ａ２  Ｔ０      ＝ｂ２  Ｓ０      …（２０−２）
…
Ｔｐ  ＋ａ１  Ｔｐ−１  ＋ａ２  Ｔｐ−２  ＋…＋ａｐ−１  Ｔ１  ＋ａｐ  Ｔ０  ＝ｂｐ  Ｓ０      …（２０−ｐ）
となる。同様にして、（１０−０）〜（１０−ｐ）式は、
Ｒ０  ＋ａ１  Ｒ１  ＋…＋ａｐ  Ｒｐ  ＝Ｔ０  ＋ｂ１  Ｔ１  ＋ｂ２  Ｔ２  ＋…＋ｂｐ  Ｔｐ      …（２１−０）
Ｒ１  ＋ａ１  Ｒ０  ＋…＋ａｐ  Ｒｐ−１  ＝ｂ１  Ｔ０  ＋ｂ２  Ｔ１  ＋…＋ｂｐ  Ｔｐ−１      …（２１−１）
Ｒ２  ＋ａ１  Ｒ１  ＋…＋ａｐ  Ｒｐ−２  ＝ｂ２  Ｔ０  ＋…＋ｂｐ  Ｔｐ−２      …（２１−２）
…
Ｒｐ  ＋ａ１  Ｒｐ−１  ＋…＋ａｐ  Ｒ０  ＝ｂｐ  Ｔ０      …（２１−ｐ）
となる。（２０−１）〜（２０−ｐ）（２１−１）〜（２１−ｐ）式に含まれるパラメータの値を推定値とし、各パラメータの推定値から構成される各式の右辺および左辺の間の残差をεｉ  とすると、
Ｔ１  ＋ａ１  Ｔ０  −ｂ１  Ｓ０  ＝ε１
Ｔ２  ＋ａ１  Ｔ１  ＋ａ２  Ｔ０  −ｂ２  Ｓ０  ＝ε２
…
Ｔｐ  ＋ａ１  Ｔｐ−１  ＋ａ２  Ｔｐ−２  ＋…＋ａｐ  Ｔ０  −ｂｐ  Ｓ０  ＝εｐ
ｂ１  Ｔ０  ＋ｂ２  Ｔ１  ＋…＋ｂｐ  Ｔｐ−１  −（Ｒ１  ＋ａ１  Ｒ０  ＋…＋ａｐ  Ｒｐ−１  ）＝εｐ＋１
ｂ２  Ｔ０  ＋…＋ｂｐ  Ｔｐ−２  −（Ｒ２  ＋ａ１  Ｒ１  ＋…＋ａｐ  Ｒｐ−２  ）＝εｐ＋２
…
ｂｐ  Ｔ０  −（Ｒｐ  ＋ａ１  Ｒｐ−１  ＋…＋ａｐ  Ｒ０  ）＝εｐ＋ｐ
が得られる。この式を更に変形すると、
【００５５】
【数２７】
2004287734.tif 000022
【００５６】
が得られる。θに対して、
【００５７】
【数２８】
2004287734.tif 000023
【００５８】
とすると、（２２）式は、（１５）式の形になる。このとき、φは、（２３）式を（６）式に適用して、
【００５９】
【数２９】
2004287734.tif 000024
【００６０】
となる。但し、Ｉはｐ×ｐの単位行列であり、
【００６１】
【数３０】
2004287734.tif 000025
【００６２】
である。（１８）式からθ（ｈ＋１）  を計算する際に、φとして（２４）式を用いることにより、（１７）式を用いる場合に比べて、（１８）式の計算が簡単になる。
【００６３】
第１発明においては、時系列をＡＲＭＡモデルに当てはめたときの前記ＡＲＭＡモデルのパラメータをガウス−ニュートン法を用いて推定する際に、白色雑音の自己相関関数、及び時系列と白色雑音との相互相関関数に理論的な値を用いることにより、確実にパラメータの値を推定することができる。
【００６４】
第２、第４、第６及び第８発明においては、白色雑音｛ｅｎ  ｝を乱数などを用いて生成し、生成した白色雑音と時系列とから、（２４）式の行列φを計算し、ガウス−ニュートン法のアルゴリズムを用いた（１８）式および（１９）式により、ＡＲＭＡモデルのパラメータの推定値を、収束するまで繰り返し計算する。
【００６５】
第３、第５及び第７発明においては、従来より確実にパラメータを推定できる方法により推定したパラメータを用いてデータを予測することにより、従来ではデータの予測ができなかった時系列についてもデータを予測することができる。
【００６６】
【発明の実施の形態】
以下本発明をその実施の形態を示す図面に基づき具体的に説明する。
図１は、本発明のデータ予測装置の構成を示すブロック図である。データ予測装置１は、本発明のパラメータ推定装置の機能を兼ね備えており、コンピュータを用いて構成されている。データ予測装置１は、演算を行うＣＰＵ（演算部）１１と、演算に伴って発生する一時的な情報を記憶するＲＡＭ（記憶部）１２と、ＣＤ−ＲＯＭドライブ等の外部記憶装置１３と、ハードディスク等の内部記憶装置１４とを備えており、ＣＤ−ＲＯＭ等の本発明に係る記録媒体２から本発明に係るコンピュータプログラム２０を外部記憶装置１３にて読み取り、読み取ったコンピュータプログラム２０を内部記憶装置１４に記憶し、ＲＡＭ１２にコンピュータプログラム２０をロードし、ＣＰＵ１１は、ロードしたコンピュータプログラム２０に基づいてデータ予測装置１に必要な処理を実行する。データ予測装置１は、キーボード又はマウス等の入力装置１５と、液晶ディスプレイ又はＣＲＴディスプレイ等の出力装置１６とを備えており、データの入力を初めとするオペレータからの操作を受け付ける構成となっている。
【００６７】
なお、データ予測装置１は、図示しない通信インタフェースを備え、通信インタフェースに接続されている図示しない外部のサーバ装置から本発明に係るコンピュータプログラム２０をダウンロードし、ＣＰＵ１１にて処理を実行する形態であってもよい。
【００６８】
図２及び３は、本発明のデータ予測装置が行うパラメータ推定の処理およびデータ予測の処理を示すフローチャートである。データ予測装置１は、オペレータの操作により入力装置１５に入力されたデータを読み込む、又は内部記憶装置１４に記憶してあるデータを読み込む等して、ＡＲＭＡモデルにおけるＡＲ過程のパラメータ（ａ１  ，ａ２  ，…，ａｐ  ）及びＭＡ過程のパラメータ（ｂ１  ，ｂ２  ，…，ｂｐ  ）の数ｐと、Ｎ個のデータから構成されるデータの時系列｛ｘｎ  ：ｎ＝１，２，…，Ｎ｝とをＲＡＭ１２に読み込む（Ｓ１）。ＣＰＵ１１は、コンピュータプログラム２０をＲＡＭ１２へロードし、ロードしたコンピュータプログラム２０に従って、ＲＡＭ１２に読み込んだ｛ｘｎ  ｝を用いて、｛ｘｎ  ｝の自己相関関数Ｒｋ  を、ｋ＝０，１，…，ｐとして、
【００６９】
【数３１】
2004287734.tif 000026
【００７０】
と計算し、（１６）式で定義される行列Ｒを生成する（Ｓ２）。ＣＰＵ１１は、次に、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、乱数生成のルーチンを用いて時刻の値を種値とする乱数を発生させる等して、Ｎ個の値から構成される平均値０の白色雑音｛ｅｎ  ：ｎ＝１，２，…，Ｎ｝を生成し（Ｓ３）、生成した白色雑音｛ｅｎ  ｝の分散Ｓ０ ２  を計算する（Ｓ４）。ＣＰＵ１１は、次に、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、データの時系列｛ｘｎ  ｝及び白色雑音｛ｅｎ  ｝を用いて、｛ｘｎ  ｝と｛ｅｎ  ｝との相互相関関数Ｔｋ  を、ｋ＝０，１，…，ｐとして、
【００７１】
【数３２】
2004287734.tif 000027
【００７２】
と計算し、（２５）式で定義される行列Ｔを生成し（Ｓ５）、計算したＳ０ ２  、生成した行列Ｒ及びＴを用いて（２４）式で定義される行列φを生成し、計算したＲｋ  及びＴｋ  を用いて（１２）式で定義される行列Ｙを生成する（Ｓ６）。ＣＰＵ１１は、次に、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、θ＝［ａ１  ，ａ２  ，…，ａｐ  ，ｂ１  ，ｂ２  ，…，ｂｐ  ］Ｔ  の初期値θ（０）  を生成し（Ｓ７）、整数ｈを、ｈ＝０とする（Ｓ８）。ＣＰＵ１１は、次に、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、ｈ番目のθの暫定値であるθ（ｈ）  を用いて、（２３）式で定義される行列ｆ（θ（ｈ）  ）を計算し（Ｓ９）、（１９）式に従って行列ε（ｈ）  を計算し（Ｓ１０）、更に、ε＝ε（ｈ）  としてＪ＝εＴε  を計算する（Ｓ１１）。
【００７３】
ＣＰＵ１１は、次に、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、Ｊの値が所定の微少値δよりも小さいか否かを判定し（Ｓ１２）、Ｊの値はδ以上であって、Ｙとｆ（θ（ｈ）  ）との残差であるεの大きさが充分な小ささではないと判定された場合は（Ｓ１２：ＮＯ）、ｈの値が予め定められている所定値よりも小さいか否かを判定する（Ｓ１３）。ｈの値が所定値よりも小さかった場合は（Ｓ１３：ＹＥＳ）、計算した行列φ及びε（ｈ）  並びにθ（ｈ）  を用いて、（ｈ＋１）番目のθの暫定値であるθ（ｈ＋１）  を（１８）式に従って計算し（Ｓ１４）、ｈに１を加えてｈをインクリメントし（Ｓ１５）、処理をステップＳ９へ戻して、Ｊの値がステップＳ１２にてδより小さいと判定されるまでステップＳ９〜Ｓ１５の処理を繰り返す。ステップＳ１３にてｈの値が所定値以上であった場合は（Ｓ１３：ＮＯ）、所定の回数だけθの暫定値の計算を繰り返してもθの値が収束しなかった場合であり、ＣＰＵ１１は、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、処理をステップＳ３へ戻して、白色雑音｛ｅｎ  ｝の生成から処理をやり直す。なお、ステップＳ１３は、計算回数ｈに基づいて判定を行うのではなく、計算時間に基づいて判定を行っても良い。
【００７４】
ステップＳ１２にてＪの値はδより小さく、Ｙとｆ（θ（ｈ）  ）との残差であるεの大きさが充分な小ささであると判定された場合は（Ｓ１２：ＹＥＳ）、現在のθの暫定値θ（ｈ）  を、パラメータθの推定値であると決定し（Ｓ１６）、決定したＡＲＭＡモデルのパラメータθ、｛ｘｎ  ｝及び｛ｅｎ  ｝を（１）式に代入して、時系列｛ｘｎ  ｝の予測データを計算し（Ｓ１７）、処理を終了する。
【００７５】
以上の処理により、ＡＲＭＡモデルのパラメータθが推定され、時系列の新たなデータが予測できる。ガウス−ニュートン法を用いた従来のパラメータ推定方法では、（１７）式を用いて計算を行っており、行列φはｐ×ｐのフルマトリックスが４つで構成されている。これに対して、本発明では、行列φは、（２４）式の如く、ｐ×ｐの一つのフルマトリックス、２つの三角行列および一つの対角行列で構成されているため、従来の方法に比べて、パラメータを計算する計算がより簡単になり、パラメータの推定値がより収束しやすくなる。また、全体の計算時間を短縮することができる。従って、従来ではデータの予測ができなかった時系列について、データの予測が可能となり、株価または売上高などの経済的な時系列について本発明を適用して、従来よりも確実に株価または売上高などのデータの予測が可能となる。
【００７６】
なお、本実施の形態においては、本発明のデータ予測装置１は、本発明のパラメータ推定装置の機能を兼ね備えている形態を示したが、これに限るものではなく、パラメータ推定装置と、ＡＲＭＡモデルによりデータを予測するデータ予測装置とを夫々別のコンピュータにて構成し、互いに情報を入出力してパラメータの推定とデータの予測とを行う形態としても良い。また、本実施の形態においては、パラメータ推定を行った後にデータ予測を行う形態を示しているが、これに限るものではなく、パラメータの推定値を記憶しておき、必要時に記憶しているパラメータの推定値を用いて予測データを計算する処理を行う形態であってもよい。
【００７７】
【発明の効果】
第１発明においては、時系列をＡＲＭＡモデルに当てはめたときの前記ＡＲＭＡモデルのパラメータを推定する際に、白色雑音の自己相関関数、及び時系列と白色雑音との相互相関関数に理論的な値を用いることにより、確実にパラメータの値を推定することができる。
【００７８】
第２、第４、第６及び第８発明においては、従来に比べて単純化された行列を含む式を用いてＡＲＭＡモデルのパラメータを推定するため、求めるべきパラメータが計算時により収束しやすくなり、従来ではパラメータの推定ができなかった時系列についても、パラメータの推定が可能となる。また、パラメータ推定の計算時間を短縮する事ができる。
【００７９】
第３、第５および第７発明においては、従来より確実にパラメータを推定できる方法により推定したパラメータを用いてデータを予測することにより、従来ではデータの予測ができなかった時系列についてもデータを予測することができる。また、株価または売上高などの経済的な時系列について本発明を適用して、従来よりも確実に株価または売上高などのデータの予測が可能となる等、本発明は優れた効果を奏する。
【図面の簡単な説明】
【図１】本発明のデータ予測装置の構成を示すブロック図である。
【図２】本発明のデータ予測装置が行うパラメータ推定の処理およびデータ予測の処理を示すフローチャートである。
【図３】本発明のデータ予測装置が行うパラメータ推定の処理およびデータ予測の処理を示すフローチャートである。
【符号の説明】
１  データ予測装置（パラメータ推定装置）
１１  ＣＰＵ（演算部）
１２  ＲＡＭ（記憶部）
２  記録媒体
２０  コンピュータプログラム
【図１】
2004287734.tif 000036
【図２】
2004287734.tif 000037
【図３】
2004287734.tif 000038