(57)【要約】【課題】制御入力を受けないARMAモデルにおいて、従来必ずしも収束性のよくなかったブートストラップ法に新たな手法を導入して計算時間の短縮を図り、的確なパラメータの推定方法を提供する。【解決手段】記憶部及び演算部を備えたコンピュータを用いて、Ｎ個のデータから構成される時系列を（ｐ、ｑ）次のＡＲＭＡモデルによって同定する際、｛ｘn｝；定常的エルゴード的正規過程ｘ（ｔ）の標本時系列（ｎ＝１、２、・・、Ｎ・・）｛ｅn｝；平均値０、分散σe2  の正規性白色雑音を用い、自己回帰（ＡＲ）部分については、３次キュムラントを用いて、パラメータを推定し、一方、移動平均（ＭＡ）部分については、2006244176.tif 000055からなるアルゴリズムを用いて推定する。【特許請求の範囲】
【請求項１】
記憶部及び演算部を備えたコンピュータを用いて、Ｎ個のデータから構成される時系列を（ｐ、ｑ）次のＡＲＭＡ（自己回帰移動平均）モデルによって同定する際、
【数１】
2006244176.tif 000042
で表され、ここで、
｛ｘn｝；定常的エルゴード的正規過程ｘ（ｔ）の標本時系列（ｎ＝１、２、・・、Ｎ・・）
｛ｅn｝；平均値０、分散σe2  の正規性白色雑音、であり、
自己回帰（ＡＲ）については、３次キュムラントを用いて、パラメータを推定し、
一方、移動平均（ＭＡ）については、数１において、簡単化してｐ＝ｑとし、
【数２】
2006244176.tif 000043
と表し、ベクトルとして、これを
2006244176.tif 000044
【数３】
2006244176.tif 000045
但し、ｓ＝σe2I      ｔ＝［Ｔ1，Ｔ2，・・，Ｔp］Tであるから、
【数４】
2006244176.tif 000046
2006244176.tif 000047
からなるアルゴリズムを用いて推定する予測のためのパラメータ推定方法。
【請求項２】
請求項１において、パラメータ推定の収束性の判定には
【数５】
2006244176.tif 000048
において、Ｄ(L)＜ε
がｍ回連続して成立したときとすることを特徴とする、予測のためのパラメータ推定方法。
【請求項３】
請求項１において、アルゴリズムを
正規性白色雑音を発生させ、｛ｅn｝の初期値とするステップ１、
2006244176.tif 000049
からなるアルゴリズムを用いて推定する予測のためのパラメータ推定方法。
【請求項４】
記憶部及び演算部を備えたコンピュータを用いて、Ｎ個のデータから構成される時系列を（ｐ、ｑ）次のＡＲＭＡ（自己回帰移動平均）モデルによってパラメータを推定し、同定する際、
N個のデータから構成される時系列｛ｘn：n=１、２、・・、N｝を記憶する第１の手段と、
平均値が０であるN個の値が時系列で並んだ白色雑音｛ｅn：n＝1、２、・・、N｝を生成し、｛ｅn｝の初期値とする第２の手段と、
2006244176.tif 000050
ＡＲＭＡモデル式から３次キュムラント式に変換し、aの推定値を計算する第４の手段と、
一方、移動平均（ＭＡ）については、数１において、簡単化してｐ＝ｑとし、
【数６】
2006244176.tif 000051
2006244176.tif 000052
【数７】
2006244176.tif 000053
2006244176.tif 000054
【請求項５】
数２より、｛ｅn｝を推定する第７の手段を用いて、前記請求項４の第３、４，５，６の手段を収束するまで繰り返す第８の手段を用いることを特徴とする請求項４記載の予測のためのパラメータ推定装置。
【請求項６】
請求項４の各手段を手順とすることを備えたことを特徴とするコンピュータプログラム。
【請求項７】
請求項５の各手段を手順とすることを特徴とするコンピュータプログラム。
【請求項８】
請求項６の各手順のコンピュータプログラムを記録したことを特徴とする記録媒体。
【請求項９】
請求項７の各手順のコンピュータプログラムを記録したことを特徴とする記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、販売予測、株価予測、炉温度制御などの予測に関するシステム同定手段を提供する。
【背景技術】
【０００２】
統計的手法による、販売予測や株価予測、工場における炉温度制御に用いるパラメータを用いて、予測するには、ＡＲＭＡ（ａｕｔｏｒｅｇｕｒｅｓｓｉｖｅ  ｍｏｖｉｎｇ  ａｖｅｒａｇｅ、自己回帰移動平均）モデルによるシステム同定手法として、拡大最小二乗法、逐次最尤法、補助変数法、疑似線形回帰法などがある。
疑似線形回帰法はブートストラップ法を用いて繰り返し計算を行うが収束しない場合が多いという欠点がある。
（ｐ、ｑ）次のＡＲＭＡモデルは、
【数１】
2006244176.tif 000002
で表され、ここで
｛ｘn｝；定常的エルゴード的正規過程ｘ（ｔ）の標本時系列（ｎ＝１、２、・・、Ｎ・・）
｛ｅn｝；平均値０、分散σe2  の正規性白色雑音、である。
【数２】
2006244176.tif 000003
とおくと
Ａ（ｚ）、Ｂ（ｚ）は既約で、定常条件、可逆条件、強正実条件を満たすものとする。
尚、数１は一般的にはＡＲＭＡＸ（Ａｕｔｏｒｅｇｒｅｓｓｉｖｅ  Ｍｏｖｉｎｇ  Ａｖｅｒａｇｅ  ｗｉｔｈ  ｅｘｇｅｎｏｕｓ  ｉｎｐｕｔ）モデルであり、
【数３】
2006244176.tif 000004
において、制御入力を受けないシステムのモデルと言い換えることができる。ここで、
｛ｕn｝は入力時系列である。
数１のパラメータ｛ａi｝、｛ｂj｝の推定に際しては、右辺そのものが有色雑音となっており、通常の最小二乗法を用いてもバイアスのかかった推定値となる。推定値がバイアスを持たないようにするために拡大最小二乗法、一般化最小二乗法、逐次最尤法、補助変数法、擬似線形回帰法等がある。しかし数３におけるパラメータ｛ａi｝、｛ｂj｝、｛ｃk｝の推定に際しては、外生入力項目部分を左辺にもっていくと、｛ｘn｝、｛ｕn｝の既知データと｛ａi｝、｛ｃk｝の未知パラメータで右辺は数１のそれと同じ形となり、ＡＲＭＡモデルのシステム同定手法の拡張と見ることができる（非特許文献１参照）。
【０００３】
制御入力を受けないARMAモデルにおいて、従来必ずしも収束性のよくなかったブートストラップ法に発明者らは新たな手法を導入してその改善を図った。即ち、モデルの構造的な前提からくる、ａ  ｐｒｉｏｒｉ  ｋｎｏｗｌｅｄｇｅ（先験的知識）を活用することによって、ロバストなパラメータ推定及び収束解を得るまでトータルな計算時間短縮が図れた（非特許文献２参照）。
コンピュータプログラムへの応用例については、特許文献１に示される。
【特許文献１】特開２００４−３０２６９０号公報
【非特許文献１】「システム同定入門」（片山徹著、朝倉書店  １９９４年）
【非特許文献２】「新しい経営・経済数学」（平成１６年４月２０日発行、株式会社中央経済社）
【発明の開示】
【発明が解決しようとする課題】
【０００４】
本発明において、ＡＲパラメータ推定は３次キュムラントを用いて推定し、ＭＡパラメータは自己相関関数及び相互相関関係を用いてモデルの構造的な前提からくるａ  ｐｒｉｏｒｉ  ｋｎｏｗｌｅｄｇｅ（先験的知識）を用いることによって推定するという方法を用いる。即ち、従来収束しにくかったケースでも収束するロバストなパラメータ推定を行い、このことにより収束解を得るまでのトータル的な計算時間短縮も図られる。
【課題を解決するための手段】
【０００５】
（ｐ、q）次のＡＲＭＡモデル（自己回帰移動平均モデル）は
【数４】
2006244176.tif 000005
で表され、確率変数が正規分布に従うとき３次以上キュムラントが０となることを利用し、式１の両辺にｘn-L2をかけて平均をとると、｛ｅn｝が３次白色であるので、
【数５】
2006244176.tif 000006
となり、更に
【数６】
2006244176.tif 000007
により、３次キュムラントとなる。数５でＬを１からｐまで動かすと、３次キュムラントに対するユール・ウオーカ方程式の
【数７】
2006244176.tif 000008
を得る。
2006244176.tif 000009
【数８】
2006244176.tif 000010
で表される。
【０００６】
一方、ユール・ウオーカ方程式から、
【数９】
2006244176.tif 000011
と表記すると、
【数１０】
2006244176.tif 000012
2006244176.tif 000013
【０００７】
移動平均ＭＡパラメータの推定について
【数１１】
2006244176.tif 000014
とおくと、
2006244176.tif 000015
【数１２】
2006244176.tif 000016
また、｛ｘn]と[ en]の相互相関関数は[ en]が白色雑音であるから、数１３となる。
【数１３】
2006244176.tif 000017
また、
【数１４】
2006244176.tif 000018
である。
尚、有限個のデータで実際に計算すると、
【数１５】
2006244176.tif 000019
のように、本来は無相関で理論的には０であるはずのものが、０でない値を持つことがある、これを
【数１６】
2006244176.tif 000020
のように表すことにする。また、数１５及び数１６より、
【数１７】
2006244176.tif 000021
と表すことにする。
【０００８】
現代制御理論において、系の状態空間表現を行い、ARMA過程表現に落し込むと、q≦p
となることは周知である。以下簡単のためq＝ｐとする。
数１において、
【数１８】
2006244176.tif 000022
と表し、ベクトルを以下のように定義する。
【数１９】
2006244176.tif 000023
【数２０】
2006244176.tif 000024
数１８は、
【数２１】
2006244176.tif 000025
と表せる。ここで、
【数２２】
2006244176.tif 000026
2006244176.tif 000027
【数２３】
2006244176.tif 000028
で与えられる。ここで、N→∞とすると、数２３の右辺は、各々｛en｝の自己相関関数と、
2006244176.tif 000029
【数２４】
2006244176.tif 000030
となる。ただし
【数２５】
2006244176.tif 000031
である。数２４を書き直しと、
【数２６】
2006244176.tif 000032
となり、これより、
【数２７】
2006244176.tif 000033
が得られる。
【０００９】
パラメータ推定の全体のアルゴリズムは、式１のものが考えられる。
式１
2006244176.tif 000034
からなるアルゴリズムを用いて推定する手段を用いる。
【００１０】
パラメータ推定の収束性の判定は
【数２８】
2006244176.tif 000035
において、D (L)＜ε
がｍ回連続して成立したときとする。
2006244176.tif 000036
に求めることができ、時間の短縮にも効果大である。
従来の方法であれば、Sの全成分が埋まった形でパラメータ推定計算がなされることになるが、本発明では、構造的な前提からくる、ａ  ｐｒｉｏｒｉ  ｋｎｏｗｌｅｄｇｅ （先験的知識）を活用して、解のロバスト、また計算時間の短縮をもできる。
【００１１】
従来手法では、数２４に相当するものとしては、
【数２９】
2006244176.tif 000037
ここで、
【数３０】
2006244176.tif 000038
である。
【００１２】
更に、本発明におけるパラメータ推定の他の全体アルゴリズムとしては、式２のものもある。
式２
2006244176.tif 000039
からなるアルゴリズムを用いて推定する。
【００１３】
解の収束性について、本発明は相互相関関数を用いてモデルの構造的な前提からくる  ａ  ｐｒｉｏｒｉ  ｋｎｏｗｌｅｄｇｅ（先験的知識） を活用することによって、ロバストなブートストラップ型アルゴリズムを開発し、これを用いて、従来手法では約半数のケースしか収束しなかったのに対し、本発明では全ケースが収束した。
【００１４】
計算の見積もりに当たって、次表のものを使用した。
表１
ＡＸ＝ｂ
ＡTＡＸ＝ＡTｂ
逆行列計算には
ガウスードウーリトル法、掃き出し法
ＱＲ分解    − Ｇｒａｍ−ｓｃｈｍｉｄｔ  法
− Ｈｏｕｓｅ  ｈｏｌｄｅｒ 法
Ｃｏｎｊｕｇａｔｅ  ｇｒａｄｉｅｎｔ  ｍｅｔｈｏｄ (共役勾配法)
などの方法がある。
ｎ変数の逆行列計算において、計算数は、
ｎ回のピボット×［１ピボットにつき｛（2n-1）回の除算＋（2n-1）回の減算｝×(n-k+1)］
となる。
除算をJ、減算をG、乗算をI、加算をK、とおくと、
2006244176.tif 000040
のオーダーとなる。
次に比較した結果を示す。
n変数の場合、双方ともS0 回繰り返しで収束と仮定し、表２の結果を得た。
【００１５】
【表１】
2006244176.tif 000041
【発明の効果】
【００１６】
従来の方法では収束しなかったケースでも、確実に収束しており、非常に優れたものになっている。本発明を活用することにより、株価予測、販売予測など確実にシステム同定でき、予測精度も高く、計算時間を従来に比べてピボット操作数でみて１０分の１前後という計算量の見積もり結果が出ている。
【発明を実施するための最良の形態】
【００１７】
以下本発明を実施の形態を示す図面に基づいて詳述する。図 1 は、本発明のパラメータ推定装置の構成を示すプロック図である。データ予測装置 1 は、本発明のパラメータ推定装置の機能を兼ね備えており、例えばパーソナルコンピュータを用いて構成されている。以下では、データ予測装置及びパラメータ推定装置をパーソナルコンピュータ 1 であるものとして説明する。
【００１８】
パーソナルコンピュータ 1 は、バス１７ を介して演算部としての CPU (  Ｃｅｎｔｒａ１  Ｐｒｏｃｅｓｓｉｎｇ  Uｎｉｔ  )１１ 、キーボード、マウス等からなる入力部 １３ , 液晶ディスプレイ等の表示部１４ 、ハードディスク等の記憶部１５、演算に伴って発生する一時的な情報を記憶する RAM (  Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ )１２ , モデム等の通信部１６, ＣＤ−ＲＯＭ  ドライブ等の外部記憶装置１８等が接続されているＣＰＵ１１ は、バス１７ を介して上述したハードウエア各部と接続されており、 CPU １１ はそれらを制御すると共に、記憶部１５ に格納された制御プログラム１５P に従って、種々のソフトウエア的機能を実行する。制御プログラム１５P は、以下に述べるソフトウエア処理を C 言語等のプログラミング言語で記述されている。
【００１９】
また、 CD 一 ROM 等の本発明に係る記録媒体２から本発明に係るコンピュータプログラム２０ を外部記憶装置１８ にて読み取り、読み取ったコンピュータブログラム２０ を記憶部１５にインストールする形態であっても良い。この場合、記憶部１５にインストールされた制御プログラムは、 CPU１１ の指示により  ＲＡＭ１２ にロードされ、  ＣＰＵｌｌはロードしたコンピュータプログラム２０ に基づいて必要な処理を実行する。通信部１６ には通信網 N （以下、インターネット N という）を介して、本発明に係るコンピュータプログラム２０ を格納した WWW (  Ｗｏｒｌｄ  Ｗｌｄｅ  Ｗｅｂ ）サーバS が接続されている。本発明にあっては、 CD 一 ROM 等の記録媒体 2 によりインストールするほか、W WW サーバS からコンピュータプログラム２０をダウンロードし、記憶部１５ にインストールする形態であってもよい。
【００２０】
ＷＷＷサーバＳにはデータの予測の基礎となる時系列データが記憶されており、パーソナルコンピュータ１はＷＷＷサーバＳから配信される時系列データを利用して将来的なデータを予測する。以下では取り扱うデータとして株価のデータに適用した場合について説明するが、これに限られるものではなく、売上高の他、データは、入力部１３からオペレータにより入力されるデータ又は図示しないセンサから出力される機器の各種データ（温度データ、流量データ）等であってもよい。
【００２１】
図２は、パーソナルコンピュータ１が行うパラメータ推定の処理及びデータ予測の処理を示すフローチャートである。パーソナルコンピュータ１のＣＰＵ１１は、ＷＷＷサーバＳからデータをダウンロードする。ＣＰＵ１１はダウンロードしたＮ個のデータから構成されるデータの時系列｛Ｘn：ｎ＝１、２、・・、Ｎ｝をＲＡＭ１２に読み込む（ステップＳ２１）。ＣＰＵ１１は、コンピュータプログラム２０をＲＡＭ１２へロードする。まず、ＣＰＵ１１は数（１４)及び数（１６）式に示されるベクトルｃ及びベクトルｃsを計算する（ステップＳ２２）．具体的には、ＣＰＵ１１はＲＡＭ１２に読み込んだ時系列データ｛Ｘｎ：ｎ＝１，２、・・、Ｎ｝を数（１５）式に代入することにより計算する。
【００２２】
続いて、ＣＰＵ１１は計算したベクトルｃ及びベクトルｃsを数（１０）式に代入することにより、ＡＲパラメータであるベクトルａを計算する（ステップＳ２３）。尚、（７）式、（８）式、（９）式及び（１０）式等の計算に必要な式は記憶部１５に記憶されており、適宜ＲＡＭ１２に読み出されてＣＰＵ１１にて計算される。次に、ＣＰＵ１１は、Ｎ個の値から構成される平均値０の白色雑音｛ｅｎ：ｎ＝１，２、・・、Ｎ｝を生成する（ステップＳ２４）。このとき、乱数生成のルーチンを用い、時刻の値を初期値として乱数を発生させるなどして、白色雑音を生成する。
【００２３】
ＣＰＵ１１は、ＲＡＭ１２にロードしたコンピュータプログラム２０に従って、生成した白色雑音｛ｅn｝の分散σe2を計算し（ステップＳ２５）、データの時系列｛Ｘａ｝及び白色雑音｛ｅn｝から、数（２５）式で定義された行列ｔと、数（１６）式で定義された行列Ｔとを計算する（ステップＳ２６）。ＣＰＵ１１は、計算した行列ｔ、Ｔ及びステップＳ２３で計算したベクトルａ、ステップＳ２５で計算した白色雑音の分散σｅ２を、数（２７）式に代入して、求めるべきＭＡパラメータであるベクトルｂを計算する（ステップＳ２７）。
【００２４】
ＣＰＵ１１は、次にＲＡＭ１２にロードしたコンピュータプログラム２０に従って、計算したベクトルａ及びｂが収束したか否かを判定する（ステップＳ２８）。ＣＰＵ１１は計算したベクトルａ及びｂが収束していないと判断した場合は（ステップＳ２８でＮＯ）、ベクトルａ及びｂとデータの時系列｛Ｘｎ｝とを（１）式にて代入して、新たな白色雑音｛ｅn｝を計算し（ステップＳ２９）、処理をステップＳ２５へ戻す。そして、ＣＰＵ１１はベクトルａ及びｂが収束するまでステップＳ２５〜２９の処理を繰り返す。尚、計算に必要な式は記憶部１５に記憶されており、適宜ＲＡＭ１２に読み出されてＣＰＵ１１にて計算される。
【００２５】
ステップＳ２８において、ＣＰＵ１１は例えば、最新のベクトルａ、ｂと前回に計算したベクトルａ、ｂであるベクトルａ'、ｂ'とから、
Δ＝│ａ−ａ'│＋│ｂ−ｂ'│
を計算し、Δの値が所定値よりも小さい場合に、ベクトルａ及びｂが収束したと判定するようにすればよい。尚、所定値は予め記憶部１５に記憶されている。
【００２６】
ステップＳ２８において、ＣＰＵ１１はベクトルａ及びｂが収束していると判断した場合（ステップＳ２８でＹＥＳ）、ＣＰＵ１１は、ＲＡＭ１２にロードしたコンピュータプログラム２０にしたがって、計算したパラメータを用いて、ＡＲＭＡモデルにより時系列｛Ｘｎ｝の予測データを計算し（ステップＳ２１０）、処理を終了する。以上の処理により、ＡＲＭＡモデルのパラメータが推定され、時系列の新たなデータが予測できる。
【図面の簡単な説明】
【００２７】
【図１】本発明のパラメータ推定装置の構成を示すブロック図である。
【図２】パーソナルコンピュータが行うパラメータ推定の処理及びデータ予測の処理を示すフローチャートである。
【符号の説明】
【００２８】
１      パーソナルコンピュータ
１１    ＣＰＵ（演算部）
１２    ＲＡＭ
１３    入力部
１４    表示部
１５    記憶部
１６    通信部
１７    バス
２０    コンピュータプログラム
Ｎ      ネットワーク
Ｓ      サーバ
【図１】
2006244176.tif 000056
【図２】
2006244176.tif 000057