JP5826893B1

JP5826893B1 - 変化点予測装置、変化点予測方法、及びコンピュータプログラム

Info

Publication number: JP5826893B1
Application number: JP2014114304A
Authority: JP
Inventors: 克彦永倉
Original assignee: Yanmar Co Ltd
Current assignee: Yanmar Co Ltd
Priority date: 2014-06-02
Filing date: 2014-06-02
Publication date: 2015-12-02
Anticipated expiration: 2034-06-02
Also published as: JP2015228183A

Abstract

【課題】将来の変化点を予測することができる変化点予測装置、変化点予測方法、及びコンピュータプログラムを提供する。【解決手段】時系列的に連続する事象に係る値及び該値に関連する複数の関連事象に係る値を複数の項目毎の項目値として前記値に対応付けて時系列に記録し、記録期間よりも後の予測対象期間における前記事象の変化点を予測する変化点予測装置であって、前記複数の項目の項目値と、該項目値に対応する期間から所定期間後の前記事象に係る値との相関の強さに基づいて予測値を算出する処理を実行し、処理により得られた予測値を記録されている実績値に連続させた第１の時系列データを作成し、該第１の時系列データに対してトレンドを算出し、トレンドの勾配を算出し、算出された勾配のマハラノビスの汎距離を算出し、算出されたマハラノビスの汎距離に基づき前記予測対象期間を含む期間における変化点を検出する。【選択図】図１７

Description

本発明は、時系列的に変化する事象を分析し将来の事象を予測する方法に関し、特に、時系列的な事象に対する値の過去の系列に基づき、該系列の将来の変化点を予測することを可能とする変化点予測装置、変化点予測方法、及びコンピュータを変化点予測装置として機能させるコンピュータプログラムに関する。

経時的に変化する事象について、該事象に係る任意の時点までの時系列データを用いて以後の変化を予測することは様々な分野で有用である。例えば為替、株価の過去の推移の観測に基づいて一日後、一週間後、半年後の為替又は株価の予測を行なったり、製品の販売実績に基づいて次期の需要予測を求めたりすることは、人間の経験及び判断力に基づいて行なわれている。上述の為替、株価、製品需要などの予測は非常に困難であるが、それらの予測が人間の主観を除去した上で高精度に実現することが可能であるとすれば、生産、販売及び在庫の管理、並びに物流、製品開発等の計画に非常に有用となる。

発明者はこれまでに、時系列的に変化する種々の事象に関する予測方法を提案してきた（特許文献１等）。既に提案した予測方法により、種々の事象の予測を実現できることが示されているが、更に予測精度を高めることが求められる。発明者が提案してきた予測方法の精度を高めるためには、過去の期間における時系列的な変化傾向（トレンド）を正確に捉えるべく、予測に用いる観測情報から適切に情報を選択することが有効である。情報を選択するための基準として、時系列的に変化する事象について変化点として定義される時点を特定することが有用である。変化点は言わば事象の潮目が変わる時点であり、観測される事象が実際に変化する時点、事象を発現させる状態が変化する時点のことである。変化点の検出は、データマイニング等の各種分野でも有用とされる。

変化点は、観測系によるノイズ、又は突発的な原因による偶然の変動（所謂外れ値）の影響を受ける。ノイズ、外れ値の除去、及び、外れ値の影響を受けない変化点の検出については、統計学、機械学習、データマイニング等の種々の分野にて多くの提案がなされている。特に時系列的な変化を考慮した変化点の検出については、時系列データの性質に対応する統計モデルを用いて外れ値及び変化点を統一的に検出する方法が特許文献２に開示されている。

特許第５０６８３８２号公報特許第３８２１２２５号公報

既に観測された時系列データに対して変化点を検出することは、過去の事象における特性を把握するために重要である。しかしながら、経時的に変化する事象についての予測で必要とされることは、過去よりも将来に起こる変化点を精度よく予測することである。

特許文献２に開示されている方法を含め従来技術では、観測された時系列データにおける言わば過去の変化点の検出は可能であるが、将来の変化点の検出は困難であった。

本発明は斯かる事情に基づいてなされたものであり、将来の変化点を予測することができる変化点予測装置、変化点予測方法、及びコンピュータを変化点予測装置として機能させるコンピュータプログラムを提供することを目的とする。

本発明に係る変化点予測装置は、時系列的に連続する事象に係る値及び該値に関連する複数の関連事象に係る値を複数の項目毎の項目値として前記値に対応付けて時系列に記録する記録手段と、該記録手段から前記値又は該値に対応する前記複数の項目毎の項目値を読み出して演算を行なう演算手段とを備え、記録期間よりも後の予測対象期間における前記事象の変化点を予測する変化点予測装置であって、前記演算手段は、前記複数の項目の項目値と、該項目値に対応する期間から所定期間後の前記事象に係る値とを対応付け、各項目について対応付けられた前記項目値と前記値との間の相関の強さを示す要因効果値に基づき選択される複数の項目の項目値から、前記事象に係る値の前記予測対象期間における予測値を算出する予測値算出手段と、記録期間における前記事象に係る値と、前記予測値算出手段により算出された前記予測対象期間における予測値とを連続させた第１の時系列データを作成する手段と、作成された前記第１の時系列データに対し、該第１の時系列データのトレンドをトレンドモデルに基づき算出するトレンド算出手段と、該トレンド算出手段が算出したトレンドの勾配を算出するトレンド勾配算出手段と、算出された勾配のマハラノビスの汎距離を算出するトレンド勾配距離算出手段と、算出されたマハラノビスの汎距離に基づき前記予測対象期間を含む期間における変化点を検出する変化点検出手段とを備えることを特徴とする。

本発明に係る変化点予測装置は、前記トレンド算出手段は、前記事象に係る値は該値を発現させる状態が前記値に与える影響を表現する関数と白色雑音との和で記述される状態空間モデルに、各時点における前記事象に係る値がトレンドと白色雑音との和で記述されるトレンドモデル、及び前後する時点におけるトレンドは略等しいとするモデルを適用して状態空間モデルの解法に基づきトレンドを算出することを特徴とする。

本発明に係る変化点予測装置は、前記勾配算出手段は、前後する時点の内、前の時点におけるトレンドから後の時点におけるトレンドへの変化率を算出するようにしてあることを特徴とする。

本発明に係る変化点予測装置は、前記演算手段は、記録期間の内、前記予測対象期間の直近の期間における前記事象に係る値の経時変化傾向と同一傾向である期間を前記記録手段に記録されている期間から特定する期間特定手段と、特定された期間における前記事象に係る値に対応する前記複数の項目毎の項目値のマハラノビスの汎距離を時系列に算出して第２の時系列データを作成する手段と、作成された前記第２の時系列データに対し、該第２の時系列データのトレンドをトレンドモデルに基づき算出する項目トレンド算出手段と、該項目トレンド算出手段が算出した項目トレンドの勾配を算出する項目トレンド勾配算出手段と、算出された勾配のマハラノビスの汎距離を算出する項目トレンド勾配距離算出手段と、算出されたマハラノビスの汎距離に基づき前記予測対象期間を含む期間における前記事象に係る値に対応する複数の項目の変化点を検出する項目変化点検出手段とを備えることを特徴とする。

本発明に係る変化点予測装置は、前記変化点検出手段は、前記トレンド勾配距離算出手段により算出されたマハラノビスの汎距離を、前記記録期間中におけるトレンド勾配のマハラノビスの汎距離の最大値で除算した第１の変化点スコアを求める手段を更に備え、前記項目変化点検出手段は、前記項目トレンド勾配距離算出手段により算出されたマハラノビスの汎距離を、前記特定された期間におけるトレンド勾配のマハラノビスの汎距離の最大値で除算した第２の変化点スコアを求める手段を更に備え、前記演算手段は、前記第１の変化点スコアと第２の変化点スコアとを加算して変化点スコアを算出するスコア算出手段を更に備えることを特徴とする。

本発明に係る変化点予測方法は、時系列的に連続する事象に係る値及び該値に関連する複数の関連事象に係る値を複数の項目毎の項目値として前記値に対応付けて時系列に記録する記録手段を備えるプロセッサが、前記記録手段から前記値又は該値に対応する前記複数の項目毎の項目値を読み出して演算を行ない、記録期間よりも後の予測対象期間における前記事象の変化点を予測する変化点予測方法であって、前記プロセッサは、前記記録手段に記録されている複数の項目の項目値と、該項目値に対応する期間から所定期間後の前記事象に係る値とを対応付け、各項目について対応付けられた前記項目値と前記値との間の相関の強さを示す要因効果値に基づき選択される複数の項目の項目値から、前記事象に係る値の前記予測対象期間における予測値を算出し、記録期間における前記事象に係る値と、算出された前記予測対象期間における予測値とを連続させた第１の時系列データを作成し、作成した第１の時系列データに対し、該第１の時系列データのトレンドをトレンドモデルに基づき算出し、算出したトレンドの勾配を算出し、算出した勾配のマハラノビスの汎距離を算出し、算出したマハラノビスの汎距離に基づき前記予測対象期間を含む期間における変化点を検出することを特徴とする。

本発明に係るコンピュータプログラムは、時系列的に連続する事象に係る値及び該値に関連する複数の関連事象に係る値を複数の項目毎の項目値として前記値に対応付けて時系列に記録する記録手段を備えるコンピュータに、前記記録手段から前記値又は該値に対応する前記複数の項目毎の項目値を読み出して演算を行ない、記録期間よりも後の予測対象期間における前記事象の変化点を予測させるコンピュータプログラムであって、前記コンピュータに、前記記録手段に記録されている複数の項目の項目値と、該項目値に対応する期間から所定期間後の前記事象に係る値とを対応付けるステップと、各項目について対応付けられた前記項目値と前記値との間の相関の強さを示す要因効果値に基づき選択される複数の項目の項目値から、前記事象に係る値の前記予測対象期間における予測値を算出するステップと、記録期間における前記事象に係る値と、算出された前記予測対象期間における予測値とを連続させた第１の時系列データを作成するステップと、作成した第１の時系列データに対し、該第１の時系列データのトレンドをトレンドモデルに基づき算出するステップと、算出したトレンドの勾配を算出するステップと、算出した勾配のマハラノビスの汎距離を算出するステップと、算出したマハラノビスの汎距離に基づき前記予測対象期間を含む期間における変化点を検出するステップとを実行させることを特徴とする。

本発明では、時系列に連続する事象に係る値について、関連項目の項目値と所定時間の時間差で対応付けられた時間差モデルにて算出される要因効果値に基づき予測値が算出され、算出された予測値が記録期間の値に連続するように第１の時系列データが作成される。作成された第１の時系列データに対し、予測期間を含む期間全体におけるトレンドが算出され、トレンドの勾配の周辺との相関強度（類似性）がマハラノビスの汎距離により算出される。変化点は類似性が低いためにマハラノビスの汎距離が大きく算出されることから、算出結果により未来を含む期間における変化点が予測可能となる。

本発明では、トレンド算出において、事象に係る値を観測系列とする状態空間モデルに、各観測系列は各時点のトレンドと白色雑音との和で記述されるというトレンドモデル、及び前後する時点のトレンドは略等しいというモデルを適用して、状態空間モデルの解法に基づきトレンドが推定算出される。これにより、トレンドの勾配はトレンド変化に敏感に反応する値として算出することが可能となり、トレンドが変化していない場合は、変化の度合いが抑えられるように算出される。

本発明では、トレンドの勾配として、前後する時点におけるトレンド間において後の時点のトレンドへの変化率が算出される。時間的に後の時点のトレンドへのトレンド増の変化が正値として算出される。

本発明では、時系列に連続する事象に係る値の関連項目のマハラノビスの汎距離が算出されて単一系列の時系列データ化され、トレンドが算出される。なお、算出対象は、予測対象期間の直近の期間における前記事象に係る値の経時変化傾向と同一傾向である期間とし、これにより事象に係る値に相関が強い項目のトレンドが求められる。項目のマハラノビスの汎距離を算出することにより、関連項目の値が影響する事象の総合的なトレンドが算出される。算出された総合的なトレンドの勾配の類似性が算出され、算出結果により項目値に基づく変化点が予測可能となる。

本発明では、事象に係る値の予測値を含む第１の時系列データから導出される変化点スコアと、前記事象に係る値の複数の項目のマハラノビスの汎距離である第２の時系列データから導出される変化点スコアとが加算されて総合的な変化点スコアとして算出される。トレンドを考慮して変化点を予測して、一時的な変動による影響を低減させた上更に、信号値の背景となる“場”の変化点が加味された変化点スコアを得ることが可能となる。

本発明による場合、将来の変化点の予測が可能である。

実施の形態１における予測装置の構成を示すブロック図である。記録部に記録される信号値及び項目毎の項目値の内容例を示す説明図である。実施の形態１における予測装置の制御部による処理手順の一例を示すフローチャートである。実施の形態１における予測装置の制御部による処理手順の一例を示すフローチャートである。時間差モデルを概念的に説明する説明図である。トレンド勾配の算出概念を示すグラフである。Ａは、実施例の信号値及び信号値に対して算出される信号トレンドを示すグラフであり、Ｂは、実施例で算出される項目ＭＤ及びＭＤトレンドを示すグラフである。Ａは、算出された項目ＭＤの正規確率及び正規確率トレンドを示すグラフであり、Ｂは、再計算された項目ＭＤ及び項目ＭＤトレンドを示すグラフである。実施例で算出された予測値を示すグラフである。実施例で算出された信号トレンド及びトレンド勾配を示すグラフである。実施例で算出されたトレンド勾配のＭＤを示すグラフである。実施例で算出された変化点スコアを示すグラフである。実施の形態２における予測装置の制御部による処理手順の一例を示すフローチャートである。実施例の項目ＭＤ、ＭＤトレンド及びトレンド勾配を示すグラフである。実施例の項目ＭＤトレンドの勾配のＭＤ及び変化点スコアを示すグラフである。実施例の複数項目に基づいて算出された変化点スコアを示すグラフである。実施例にて得られる信号値に基づく変化点スコアと、項目値に基づく変化点スコアとを合わせて表示したグラフである。

以下、本発明の実施の形態を、図面に基づいて具体的に説明する。

（実施の形態１）
図１は、実施の形態１における予測装置１の構成を示すブロック図である。予測装置１は例えばパーソナルコンピュータ、サーバコンピュータ等のコンピュータを用いる。予測装置１は、制御部１０、記録部１１、一時記憶部１２、入力部１３及び出力部１４を備える。

制御部１０は、ＣＰＵ（Central Processing Unit ）を用いる。制御部１０は、以下に説明する予測プログラム１Ｐに基づき、パーソナルコンピュータを制御し、本実施の形態における予測装置１としての機能を発揮させる。

記録部１１は、ＲＯＭ（Read Only Memory）、ハードディスクドライブ等の不揮発性メモリを用いる。なお記録部１１は、外付けのハードディスクドライブ、光学ディスクドライブ、通信網を介して接続される他の記録装置であってもよい。すなわち、記録部１１とは、制御部１０からアクセス可能な１又は複数の情報記録媒体の総称である。

記録部１１には、実施の形態１の変化点予測方法を実現するための各種手順を含む予測プログラム１Ｐが記録されている。また、記録部１１の記録領域の一部は、信号値及び該値に対応する複数の項目のデータ（各項目の項目値）を記録するデータベース（ＤＢ）１１０として用いられる。制御部１０は、データベース１１０に対し、信号値及び各項目の項目値の読み書きが可能である。データベース１１０は、信号値及び各項目の項目値を時系列に記録している。

一時記憶部１２は、例えば、ＤＲＡＭ（Dynamic Random Access Memory）、ＳＲＡＭ（Static Random Access Memory）等の不揮発性メモリである。一時記憶部１２は、制御部１０の処理によって発生した情報を一時的に記憶する。

入力部１３は、キーボード、マウス等を用い、ユーザの操作に基づき情報を入力する。

出力部１４は、液晶モニタなどの表示部、又はプリンタ等の印刷部を用い、制御部１０による情報の処理結果を出力する。

図２は、記録部１１に記録される信号値及び項目毎の項目値の内容例を示す説明図である。図２において、メンバーとは、単位期間毎のデータであることを示すインデックスであり、１，２，…，l として示している。信号値自体は、信号値Ｍ₁，Ｍ₂，…，Ｍ_lとして記録される。「項目１」，「項目２」，…，「項目ｋ」は、信号値Ｍ₁，Ｍ₂，…，Ｍ_lに関連する項目であり、Ｘ₁₁，Ｘ₁₂，…，Ｘ_1kは、メンバー１の信号値Ｍ₁に関連する各項目の項目値を示している。なお、信号値Ｍ及び各項目の項目値Ｘは、実際の値を用いてもよいが、入力後又は入力前に規準化しておいてから用いることが好ましい。規準化の方法は例えば、制御部１０が信号値及び項目毎の項目値を入力した後に、各項目の項目値Ｘ（例えば月平均）から、当該項目の項目値の平均（用いる全メンバーの月平均の平均値）を減算しておくなどの方法がある。

このように構成される予測装置１の制御部１０により実行される変化点予測処理について説明する。なお、以下に説明する手順には、Ｔ法に基づく予測値算出処理が含まれるが、予測値算出処理については、特許第５０６８３８２号公報、特許第５０６８３８３号公報、特許第５４１６８０９〜５４１６８１１号公報、特願２０１３−１２５９３３〜２０１３−１２５９３６等の先行文献に記載されているため、詳細な説明を一部省略する。

図３及び図４は、実施の形態１における予測装置１の制御部１０による処理手順の一例を示すフローチャートである。

制御部１０は、入力部１３から、時系列の信号値及びこれに関連する各項目の時系列の項目値を入力し、入力した信号値及び項目毎の項目値を対応付けて記録部１１のデータベース１１０に記録する（ステップＳ１０１）。

制御部１０は、記録部１１のデータベース１１０に記録した信号値及び項目毎の項目値に基づいて、時間差モデルを生成する（ステップＳ１０２）。

ステップＳ１０２において制御部１０は、時系列に記録された信号値と、項目毎の項目値とを所定期間分ずらして対応付ける。時間差モデルとは、項目毎の項目値と該項目の所定期間後の信号値とを対応付けたモデルである。図５は、時間差モデルを概念的に説明する説明図である。図５では、図の左から右へ向かって時間の経過を示している。図４の下部の各矩形は、１年毎の各項目の項目値を示し、上部の各矩形は、各項目と同時期の１年毎の信号を示している。時間差モデルでは、同時期の信号と各項目の項目値とを対応付けるのではなく、図２で示した信号値Ｍ_iに対し、所定期間前のメンバー（i-t ）の各項目の項目値Ｘ_{i-t 1}，Ｘ_{i-t 2}，…，Ｘ_{i-t k}を対応付ける。図５の例では、所定期間とは１年であり、２０００年１月の信号値Ｍに対し、１年前の１９９９年１月の項目値Ｘ₁，Ｘ₂，…，Ｘ_kを対応付け、２００７年１２月の信号値Ｍに対しては１年前の２００６年１２月の項目値Ｘ₁，Ｘ₂，…，Ｘ_kを対応付けた時間差モデルを生成する。

制御部１０は、生成した時間差モデルにおいて対応する信号値と項目毎の項目値との関係に基づき、信号値及び項目毎の項目値を変数変換する（ステップＳ１０３）。ステップＳ１０３における変数変換は、実施の形態１では、信号値について対数変換を行うと共に、対数変換した信号値に対して項目値を２次変数変換する。

制御部１０は、変数変換後の信号値及び項目毎の項目値を用いて、予測に用いる信号期間を選択する処理を行なうべく、対数変換後の信号値を用いて信号トレンドを算出する（ステップＳ１０４）。信号トレンドについては、状態空間モデルによる２次のトレンド成分を用いた計算手法を採用して算出する。

制御部１０は、算出した信号トレンドに基づき、予測対象期間の直近の信号トレンドと同じ信号区間を選択する（ステップＳ１０５）。ここで制御部１０は、予測対象期間の直近の信号トレンドを参照し、予測対象期間の直近の信号トレンドが増加傾向（＋）である場合、記録期間中において増加傾向（＋）の区間を選択し、直近の信号トレンドが減少傾向（−）である場合、記録期間において減少傾向（−）の区間を選択する。

制御部１０は、ステップＳ１０５で選択した信号区間を連続区間として用い、この期間における項目のマハラノビスの汎距離（以下、ＭＤ：Maharanobis’s Distance と言う）及びＭＤトレンドを算出する（ステップＳ１０６）。ステップＳ１０６においては、予測対象期間における信号値（未知信号）に対応付けるべき項目値をも含めて項目のＭＤ（以下、項目ＭＤと言う）を算出する。

制御部１０は、ステップＳ１０６で算出した項目ＭＤの正規確率及び該正規確率のトレンドを算出する（ステップＳ１０７）。制御部１０は、算出した正規確率及び該正規確率のトレンドに基づき、予測対象区間の直近の項目ＭＤを含む１循環分の信号区間を選択する（ステップＳ１０８）。

ステップＳ１０８において制御部１０は、以下のように信号区間を選択する。ステップＳ１０７で算出した項目ＭＤの正規確率は０〜１の間の数値であるから、制御部１０は、閾値を０．５に設定して正規確率が０．５以上となる第１トレンド区間と、０．５未満の第２トレンド区間とを判別する。制御部１０は、予測対象期間の最新の区間が第１トレンド区間であるか、第２トレンド区間であるかを判断し、第１トレンド区間である場合には直近の第１トレンド区間と１つ前の第２トレンド区間とを１循環分の信号区間として選択する。制御部１０は、予測対象区間の最新の区間が第２トレンド区間である場合には、直近の第２トレンド区間と１つ前の第１トレンド区間とを１循環分の信号区間として選択する。

次に制御部１０は、ステップＳ１０８で選択した信号区間での未知項目（予測対象期間に信号値に対応する項目値）を含む項目ＭＤ及び該項目ＭＤのトレンドを再計算する（ステップＳ１０９）。制御部１０は、未知項目ＭＤと同一符号の区間を選択する（ステップＳ１１０）。ステップＳ１１０において制御部１０は、予測対象期間の直近においてステップＳ１０９で再計算したトレンド曲線が示すＭＤと閾値（例えば１．０）とを比較し、直近のトレンド曲線が示すＭＤが閾値以上である場合、閾値以上の区間を選択し、直近のトレンド曲線が示すＭＤが閾値未満である場合、閾値未満の区間を選択する。ここで比較した閾値には全ＭＤの平均である１．０の値が用いられる。

制御部１０は、ステップＳ１１０で選択した区間における信号値と、時間差モデルにて対応する各項目の項目値との関係を特定する特定式をＴ法に基づいて算出する（ステップＳ１１１）。ステップＳ１１１において具体的には、制御部１０は信号値に対して項目毎の要因効果の強さを示す比例定数β及びＳＮ比η（２乗比）を下記の式（１）及び式（２）を適用して算出する。

制御部１０は、ステップＳ１１１にて算出した項目毎の比例定数βを用い、ステップＳ１１０で選択した信号期間における各信号値について、以下に示す式（３）により出力の推定値を算出する（ステップＳ１１２）。

次に制御部１０は、推定値についての推定精度であるＳＮ比ηを重みづけ係数として用いた総合推定値を以下に示す式（４）に基づいて算出する（ステップＳ１１３）。

制御部１０は、信号値及び総合推定値に基づいて、各項目の総合推定値のＳＮ比（ｄｂ）を以下に示す式（５）により算出する（ステップＳ１１４）。

制御部１０は、各項目について要因効果値を算出する（ステップＳ１１５）。要因効果値は、各項目について当該項目を除いた各項目の項目値の総合推定値のＳＮ比と、当該項目を含む各項目の項目値の総合推定値のＳＮ比の差分を求めることによって算出される。より具体的には、制御部１０は、各項目について当該項目を含む各項目の項目値の総合推定値のＳＮ比から、当該項目を含まない各項目の項目値の総合推定値のＳＮ比を減じた値を算出する。

制御部１０は、算出した要因効果値が大きい順に項目を複数選択し、選択した複数の項目の項目値に対する総合推定値のＳＮ比を、項目数毎に算出する（ステップＳ１１６）。制御部１０は、項目数毎の総合推定値のＳＮ比に基づき、当該ＳＮ比を最大にする項目数を決定する（ステップＳ１１７）。

制御部１０は、ステップＳ１１７にて決定した項目数分、要因効果値が大きいものから項目を選択する（ステップＳ１１８）。制御部１０は、時間差モデルにおいて予測対象期間に対応する項目の内、ステップＳ１１８で選択した項目の項目毎の項目値を式（４）に当てはめて予測値を算出する（ステップＳ１１９）。ここで制御部１０は、式（４）における比例定数β及びＳＮ比ηとして、ステップＳ１１１で算出したものを用いる。ステップＳ１１９において制御部１０は、信号値が規準化されている場合には逆変換してもよい。

ステップＳ１１９で算出された予測対象期間における予測値を、記録期間における既知の信号値と合わせて時系列データを作成する（ステップＳ１２０）。

制御部１０は、作成した時系列データに対し、隣り合う時点のトレンド成分は略等しいというトレンドモデルに基づいて予測対象期間を含む期間全体における信号トレンドを算出する（ステップＳ１２１）。ステップＳ１２１における信号トレンドの算出には、状態空間モデルによる１次のトレンド成分を適用する。具体的には制御部１０は、トレンドｔ_nについて、以下に示す式（６）にて表現されるようにトレンドが時間軸において局所的に略一定の値をとり、前のタイミングにおけるトレンドｔ_n-1に対してｔ_n≒ｔ_n-1であるというモデル（ランダムウォークモデル、参考：トレンド成分モデル−構造の確率的変化のモデル−、北川源四郎、１９９３、「ＦＯＲＴＲＡＮ７７時系列解析プログラミング」、pp.252-263）を、ステップＳ１２０で作成した時系列データを観測値とする状態空間モデルに当てはめ、カルマンフィルタ及び平滑化処理等の処理を実行することによってトレンドを算出する。

制御部１０は、ステップＳ１２１で算出したトレンドに対し、風上差分からのトレンド勾配を算出する（ステップＳ１２２）。ステップＳ１２２において詳細には、制御部１０は、各時点について算出されたトレンドに対し、対象としている一の時点の次の時点におけるトレンドへの変化率を算出する。図６は、トレンド勾配の算出概念を示すグラフである。横軸は時間軸を示し、縦軸は算出されたトレンドを示す。ステップＳ１２２にて算出されるトレンド勾配は、図６中のＧに対応する。

制御部１０は、ステップＳ１２２で算出したトレンド勾配についてＭＤを算出する（ステップＳ１２３）。これにより、トレンド勾配の類似性が求められる。なお、算出されたＭＤが近いほどトレンドが類似し、遠いほどトレンドが非類似であることが示される。

制御部１０は、ステップＳ１２３で算出したＭＤに基づき、変化点スコアを算出する（ステップＳ１２４）。ステップＳ１２４において制御部１０は、ステップＳ１２０で作成した時系列データの予測対象期間を含む期間全体に亘って算出したＭＤの最大値で、算出されたＭＤ全てを除算する。ステップＳ１２４の処理により、予測対象期間を含む全体で算出されたＭＤの内の最大値に対応する時点の変化点スコアが１．０となるように、各時点での変化点スコアが算出される。これにより、非類似性が高いトレンドの勾配に係る時点が変化点として検出される。

このようにして、予測された予測値を含む全期間でのトレンドに基づき、予測対象期間を含む期間における信号値の変化点を検出することができる。検出装置１により導出される変化点スコアは、予測対象期間を含む期間のトレンドが考慮されているから、トレンドに従って真に状態が変化することに応じて導出される。

なお、図３及び図４のフローチャートに示した処理手順における予測値算出処理（Ｓ１０２〜Ｓ１１９）の手順は一例であって本発明に係る発明者が提案する時間差モデルに基づく他の予測方法（先行技術文献特許第５０６８３８２号公報、特許第５０６８３８３号公報、特許第５４１６８０９〜５４１６８１１号公報、特願２０１３−１２５９３３〜２０１３−１２５９３６等）のいずれか及び他の公知の方法を含む方法を用いるようにしてもよい。

次に、実施の形態１の予測装置１による変化点予測方法を適用した具体的な実施例について説明する。

実施例として、実施の形態１の変化点予測方法を、為替レートに適用した例を説明する。実施例では、例えば１５の各国通貨の対ドル為替レートの２０００年１月〜２００７年１２月までの各月の平均を項目として用い、その内のポンドドル（ＧＢＰ／ＵＳＤ）を信号値として変化点予測方法に適用する。実施例における時間差モデルは、図５に示したように所定時間を１年として生成される。

予測装置１の制御部１０は、時間差モデルを生成し（Ｓ１０２）、変数変換を行なった後（Ｓ１０３）、変換後の信号値を用いて信号トレンドを算出する（Ｓ１０４）。

図７（Ａ）は、実施例の信号値及び信号値に対して算出される信号トレンドを示すグラフである。図７（Ａ）の横軸は年月を時系列に示し、縦軸は信号値の大きさを示している。なお、図７（Ａ）中では、信号値を白丸で示し、信号値に対して算出される信号トレンドを実線で示している。また図７（Ａ）には、信号トレンドの符号の信号区間を（＋）と（−）とで区別して示している。制御部１０は、直近の信号トレンドは増加傾向（＋）であるから、制御部１０は図７（Ａ）中の（＋）で示される増加傾向の信号区間を選択する（Ｓ１０５）。

次に予測装置１の制御部１０は、選択した信号区間を連続区間のように扱い、当該区間の信号値に対応する項目毎の項目値のＭＤ及びＭＤトレンドを算出する（Ｓ１０６）。

図７（Ｂ）は、実施例で算出される項目ＭＤ及びＭＤトレンドを示すグラフである。図７（Ｂ）の横軸は年月を時系列に示し、縦軸は項目ＭＤの大きさを示している。図７（Ｂ）では項目ＭＤは丸印で、ＭＤトレンドは実線にて示されている。なおステップＳ１０６では、予測対象期間における信号値（未知信号）に対応付けるべき項目値をも含めて項目ＭＤが算出されている。予測対象期間に対応する項目ＭＤは、図７（Ｂ）中で白丸により示している。

制御部１０は、項目ＭＤの正規確率及び正規確率のトレンドを算出し（Ｓ１０７）、算出された正規確率のトレンドに基づいて、直近項目ＭＤを含む１循環分の信号区間を選択する（Ｓ１０８）。

図８（Ａ）は算出された項目ＭＤの正規確率及び正規確率トレンドを示すグラフである。図８（Ａ）の横軸は年月を時系列に示し、縦軸は正規確率を示している。図８（Ａ）では、項目ＭＤの正規確率を丸印で示し、正規確率のトレンドを細い実線で示し、太い実線にて正規確率のトレンドが閾値０．５以上である第１トレンド区間と閾値０．５未満である第２トレンド区間とを区別して示している。制御部１０は、予想対象期間の最新の区間が第２トレンド区間であると判断し、直近の第２トレンド区間と１つ前の第１トレンド区間とを１循環分の信号区間として選択する。図８（Ａ）中、選択された１循環分の信号区間を「ＭＤ１循環」として示している。

制御部１０は、ステップＳ１０８で選択した信号区間での未知項目を含む項目ＭＤ及び該項目ＭＤのトレンドを再計算し（Ｓ１０９）、未知項目ＭＤと同一符合の既知信号を選択する（Ｓ１１０）。

図８（Ｂ）は、再計算された項目ＭＤ及び項目ＭＤトレンドを示すグラフである。図８（Ｂ）の横軸は年月を時系列に示し、縦軸は項目ＭＤの大きさを示している。図８（Ｂ）中、各項目ＭＤを丸印で示し、ＭＤトレンドを細い実線で示し、太い実線にてトレンドの符号が同一の区間を（＋）と（−）とで区別して示している。予測対象期間に対応する項目ＭＤは、図８（Ｂ）中で白丸により示している。制御部１０は、未知項目ＭＤの符号を（−）と判断し、（−）の符号を示すトレンドに対応する区間を選択する。図８（Ｂ）中、ステップＳ１１０で選択される区間を「未知ＭＤトレンド選択」として示しており、具体的には２００７年１月から２００７年１２月が選択されている。

制御部１０は、ステップＳ１１０にて図８（Ｂ）で示したように選択した区間（２００７年１月から２００７年１２月）における信号値及び該信号値と時間差モデルで対応する各項目の項目値との関係と特定する特定式、即ち比例定数β及びＳＮ比ηを算出する（Ｓ１１１）。制御部１０は、Ｓ１１０で選択した区間における各信号値について、算出された比例定数β及びＳＮ比ηを用いて出力の推定値を算出し、要因効果値、総合推定値のＳＮ比に基づいて予測値を算出する（Ｓ１１２〜Ｓ１１９）。

図９は、実施例で算出された予測値を示すグラフである。図８の横軸は、ステップＳ２０９で選択された区間（２００７年１月から２００７年１２月）以後の予測対象期間（２００８年１月〜１２月）を含む時間軸を示している。図９の縦軸は、ＵＳＤ／ＧＢＰのレートを示している。図９中、二重線が予測値を結ぶ線である。なお図９中、実績値を丸印で示し、白丸にて予測対象期間を含む以後の実績値を比較のために示している。

次に制御部１０は、予測値算出処理によって得られた予測対象期間における予測値（図９の二重線）を、記録期間における既知の信号値と合わせて時系列データを作成する（Ｓ１２０）。ステップＳ１２０において制御部１０は、記録期間（２０００年１月〜２００７年１２月）における信号値に、予測対象期間（２００８年１月〜１２月）の予測値を連続させて時系列データを作成する。

制御部１０は、作成された時系列データに対し、予測対象期間を含む期間全体における信号トレンドを算出し（Ｓ１２１）、トレンド勾配を算出する（Ｓ１２２）。

図１０は、実施例で算出された信号トレンド及びトレンド勾配を示すグラフである。図１０の横軸は年月を時系列に示し、左側の縦軸は信号トレンドに対応するＵＳＤ／ＧＢＰのレートを示し、右側の縦軸はトレンド勾配（信号トレンドの差分）に対応する数値を示している。なお図１０では、図示を明確とするために記録期間中、２０００年１月から２００４年１２月までの期間の信号トレンド及びトレンド勾配を省略して示している。

図１０では、変化点において値が増加する場合は勾配が正、値が減少する場合は勾配が負で示されていることがわかる。

次に制御部１０は、トレンド勾配のＭＤを算出する（Ｓ１２３）。

図１１は、実施例で算出されたトレンド勾配のＭＤを示すグラフである。図１１の横軸は図１０と対応する年月を時系列に示し、縦軸（左側）はＭＤを示す。なお図１１には比較のために図１０の信号トレンドの勾配を下部に示し、縦軸（右側）にトレンド勾配に対応する値を示している。図１１から、２００８年後半に、ＭＤが大きい値を示すことが読み取れる。

制御部１０は、トレンド勾配のＭＤを最大値で除算して変化点スコアを算出する（Ｓ１２４）。

図１２は、実施例で算出された変化点スコアを示すグラフである。図１２の横軸は図１１と対応する年月を時系列に示し、縦軸（右側）は変化点スコアを示す。図１２では、変化点スコアは白抜きの棒グラフにより示されている。図１２には比較のために信号値及び予測値を上部に示し、縦軸（左側）に信号値即ちＵＳＤ／ＧＢＰのレートを示している。図１２中、二重線は信号値及び予測値からなる時系列データに対して求められる信号値トレンド（Ｓ１２１）を示している。このようにして、予測対象期間を含む期間の時系列データに対してトレンドを考慮した変化点が検出され、２００８年後半に、ＵＳＤ／ＧＢＰのレートの変化点が予測されている。

（実施の形態２）
実施の形態１では、信号値の変化点を信号値の単一系列に基づいて予測する構成とした。これに対し実施の形態２では、信号値の予測に用いた複数の項目の総合的な変化点を合わせて予測し、信号値単独及び複数の項目の両面から信号値の変化点を予測する。

実施の形態２における予測装置１のハードウェア構成は、実施の形態１における予測装置１と同一であるから、同一の符号を付して詳細な説明を省略し、以下に、実施の形態２における予測装置１の制御部１０により実行される処理手順について詳細を説明する。

図１３は、実施の形態２における予測装置１の制御部１０による処理手順の一例を示すフローチャートである。なお図１３のフローチャートに示す処理手順の内、実施の形態１における図３及び図４のフローチャートと共通する手順については、同一のステップ番号を付して詳細な説明を省略する。

実施の形態２では、ステップＳ１０３における変数変換は、以下のように行う。変数変換は、生成した時間差モデルにおいて対応する信号値と項目毎の項目値との関係に基づき、信号値及び項目毎の項目値を変数変換するものである。しかしながら、ステップＳ１０２で生成した時間差モデルにおいて、未知信号に対応する未知項目には、対応する信号値がない。実施の形態２におけるステップＳ１０３では、制御部１０は既知信号と既知項目との関係から定められる係数を項目毎に保存しておき、未知項目の各項目に乗ずることで変数変換を行う。

制御部１０は、ステップＳ１０５で選択された区間を連続区間として、項目ＭＤを未知信号に対応する項目も含めて算出し（ステップＳ１３１）、ステップＳ１０５で選択された信号区間を連続区間とし、更に予測対象期間を接続させた期間における時系列のＭＤを時系列データとして作成し（ステップＳ１３２）、作成したＭＤの時系列データに対し、隣り合う時点のＭＤのトレンド成分は略等しいというトレンドモデルに基づいて前記期間全体におけるＭＤトレンドを算出する（ステップＳ１３３）。

ステップＳ１３３におけるＭＤトレンドの算出には、状態空間モデルによる１次のトレンド成分を適用する。具体的には制御部１０は、ＭＤトレンドＭＤＴ_nについて、式（７）にて表現されるようにＭＤトレンドが時間軸において局所的に略一定の値をとり、前のタイミングにおけるトレンドＭＤＴ_n-1に対してＭＤＴ_n≒ＭＤＴ_n-1であるというモデル（ランダムウォークモデル）を、ステップＳ１３１で作成された時系列のＭＤを観測値とする状態空間モデルに当てはめ、カルマンフィルタ及び平滑化処理等の処理を実行することによってＭＤトレンドを算出する。

次に制御部１０は、ステップＳ１３３で算出したＭＤトレンドに対して風上差分からのトレンドの勾配を算出する（ステップＳ１３４）。制御部１０は、算出した勾配のＭＤを算出し（ステップＳ１３５）、算出したＭＤを、ＭＤの最大値で除算して変化点スコアを算出し（ステップＳ１３６）、処理を終了する。

実施の形態２の予測装置１による変化点予測方法を適用した具体的な実施例について説明する。実施例として、実施の形態１同様に変化点予測方法を為替レートに適用した例を説明する。

予測装置１の制御部１０は、項目ＭＤを未知信号に対応する項目も含めて算出し（Ｓ１３１）、予測対象期間を含む期間における時系列のＭＤを時系列データとして作成し（Ｓ１３２）、ＭＤトレンドを算出し（Ｓ１３３）、ＭＤトレンドの勾配を算出する（Ｓ１３４）。

図１４は、実施例の項目ＭＤ、ＭＤトレンド及びトレンド勾配を示すグラフである。図１４の横軸は年月を時系列に示し、縦軸（左側）は項目ＭＤに対応する数値を示している。図１４中では、項目ＭＤを黒丸で示し、予測対象期間に対応する項目ＭＤを白丸で示し、ＭＤトレンドを二重線にて示している。また、図１４には下部に、算出されたトレンド勾配を示し、右側の縦軸にトレンド勾配に対応する数値が示されている。なお、図１４中の項目ＭＤ、ＭＤトレンド及びトレンド勾配は、ステップＳ１０５で選択された信号区間（図７（Ａ）参照）に対応する項目について算出されているため途切れて示されている。図１４に示すように、未知信号に対応する未知項目についても連続区間として扱われてＭＤトレンドが算出されている。ＭＤトレンドの勾配の算出により、直近の信号トレンドと類似する区間にて該区間の信号値を発現させる“場”のトレンドの局所的な変化が得られる。

予測装置１の制御部１０は、算出したＭＤトレンドの勾配のＭＤを算出し（Ｓ１３５）、算出したＭＤに基づき変化点スコアを算出する（Ｓ１３６）。

図１５は、実施例の項目ＭＤトレンドの勾配のＭＤ及び変化点スコアを示すグラフである。図１５の横軸は年月を時系列に示し、縦軸（左側）はＭＤに対応する数値を示し、縦軸（右側）は変化点スコアに対応する数値を示している。ＭＤトレンドの勾配のＭＤを算出することにより、“場”のトレンドの局所的な変化の類似性を示す指標（値が大きいほど遠く、非類似を示す）が得られる。そして変化点スコアはＭＤの最大値で除算されて算出されているから、変化点スコアの最大値は１．０である。

図１６は、実施例の複数項目に基づいて算出された変化点スコアを示すグラフである。図１６の横軸は年月を時系列に示し、縦軸（右側）に変化点スコアを示している。図１６中、変化点スコアは黒塗りの棒グラフにより示されている。図１６には比較のために図１４の項目ＭＤ及びＭＤトレンドを上部に示し、縦軸（左側）に項目ＭＤに対応する値を示している。図１６に示すように、信号値に関連する複数の項目の項目値に基づいて予測される変化点も２００８年後半に検出されている。

実施の形態１及び２に示した処理手順を組み合わせ、変化点を評価することが可能である。具体的には、図３及び図４のフローチャートに示した処理手順の内のステップＳ１０５までの処理を行なった時点で、図１３のフローチャートに示したステップＳ１３１〜Ｓ１３６を実行し、その後又は同時並行的に、図３及び図４のステップＳ１０６以降の処理を別途行なう。これにより、信号値の単一系列に基づいて予測される変化点と、信号値に関連する複数の項目の項目値に基づいて総合的に予測される変化点とを合わせて取得することができる。

図１７は、実施例にて得られる信号値に基づく変化点スコアと、項目値に基づく変化点スコアとを合わせて表示したグラフである。図１７の横軸は年月を時系列に示し、縦軸（左側）に信号値即ちＵＳＤ／ＧＢＰのレートを示し、縦軸（右側）に変化点スコアに対応する数値を示している。図１７では、白丸で予測対象期間（２００８年１月〜１２月）中の予測値を示している。予測対象期間中も含め、黒丸は実績値を示している。また図１７では、予測装置１により、信号値に基づいて算出された変化点スコアと信号値に関連する複数の項目の項目値に基づいて算出された変化点スコアとの和を白塗りの棒グラフで示し、項目に基づく変化点スコアを黒塗りの棒グラフで示している。

図１７に示すように、信号値に基づく変化点スコアと、項目に基づく変化点スコアとに基づいて、２００８年後半において変化点スコアの和「２」が算出されており、過去最大の変化点が予測されている。２００８年後半には、所謂リーマンショックによって実績値に示すようにＵＳＤ／ＧＢＰレートには大きな変化があったのであるが、本発明に係る予測装置１の処理により、対応する変化点が予測されていることがわかる。詳細に分析すれば、信号値に基づく変化点スコアは、２００８年４月から僅かに大きく予測されているが、項目に基づく変化点スコアは、２００７年７月付近から大きく検出されている。２００８年５月からは信号値に基づく変化点スコア及び項目に基づく変化点スコアのいずれにおいても変化点スコアが増幅し、２００８年９月に最大値を迎えることが予測されている。

このようにして、信号値のみならず、該信号値を発現させる関連項目の項目値から推測される“場”の変化点を、将来を含む期間全体のトレンドに基づいて予測することができる。トレンドを考慮して変化点を予測することができるから、一時的な変動による影響を低減することができ、更に、信号値の背景となる“場”の変化点を関連項目の項目値に基づいて予測することができる。これにより、信号値の将来の変化点を精度よく予測することが可能となる。

なお、上述のように開示された本実施の形態はすべての点で例示であって、制限的なものではないと考えられるべきである。本発明の範囲は、上記した意味ではなく、特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。

１予測装置
１０制御部
１１記録部
１３入力部
１Ｐ予測プログラム

Claims

時系列的に連続する事象に係る値及び該値に関連する複数の関連事象に係る値を複数の項目毎の項目値として前記値に対応付けて時系列に記録する記録手段と、該記録手段から前記値又は該値に対応する前記複数の項目毎の項目値を読み出して演算を行なう演算手段とを備え、記録期間よりも後の予測対象期間における前記事象の変化点を予測する変化点予測装置であって、
前記演算手段は、
前記複数の項目の項目値と、該項目値に対応する期間から所定期間後の前記事象に係る値とを対応付け、各項目について対応付けられた前記項目値と前記値との間の相関の強さを示す要因効果値に基づき選択される複数の項目の項目値から、前記事象に係る値の前記予測対象期間における予測値を算出する予測値算出手段と、
記録期間における前記事象に係る値と、前記予測値算出手段により算出された前記予測対象期間における予測値とを連続させた第１の時系列データを作成する手段と、
作成された前記第１の時系列データに対し、該第１の時系列データのトレンドをトレンドモデルに基づき算出するトレンド算出手段と、
該トレンド算出手段が算出したトレンドの勾配を算出するトレンド勾配算出手段と、
算出された勾配のマハラノビスの汎距離を算出するトレンド勾配距離算出手段と、
算出されたマハラノビスの汎距離に基づき前記予測対象期間を含む期間における変化点を検出する変化点検出手段と
を備えることを特徴とする変化点予測装置。
前記トレンド算出手段は、
前記事象に係る値は該値を発現させる状態が前記値に与える影響を表現する関数と白色雑音との和で記述される状態空間モデルに、各時点における前記事象に係る値がトレンドと白色雑音との和で記述されるトレンドモデル、及び前後する時点におけるトレンドは略等しいとするモデルを適用して状態空間モデルの解法に基づきトレンドを算出する
ことを特徴とする請求項１に記載の変化点予測装置。
前記勾配算出手段は、
前後する時点の内、前の時点におけるトレンドから後の時点におけるトレンドへの変化率を算出するようにしてある
ことを特徴とする請求項１又は２に記載の変化点予測装置。
前記演算手段は、
記録期間の内、前記予測対象期間の直近の期間における前記事象に係る値の経時変化傾向と同一傾向である期間を前記記録手段に記録されている期間から特定する期間特定手段と、
特定された期間における前記事象に係る値に対応する前記複数の項目毎の項目値のマハラノビスの汎距離を時系列に算出して第２の時系列データを作成する手段と、
作成された前記第２の時系列データに対し、該第２の時系列データのトレンドをトレンドモデルに基づき算出する項目トレンド算出手段と、
該項目トレンド算出手段が算出した項目トレンドの勾配を算出する項目トレンド勾配算出手段と、
算出された勾配のマハラノビスの汎距離を算出する項目トレンド勾配距離算出手段と、
算出されたマハラノビスの汎距離に基づき前記予測対象期間を含む期間における前記事象に係る値に対応する複数の項目の変化点を検出する項目変化点検出手段と
を備えることを特徴とする請求項１乃至３のいずれか１つに記載の変化点予測装置。
前記変化点検出手段は、前記トレンド勾配距離算出手段により算出されたマハラノビスの汎距離を、前記記録期間中におけるトレンド勾配のマハラノビスの汎距離の最大値で除算した第１の変化点スコアを求める手段を更に備え、
前記項目変化点検出手段は、前記項目トレンド勾配距離算出手段により算出されたマハラノビスの汎距離を、前記特定された期間におけるトレンド勾配のマハラノビスの汎距離の最大値で除算した第２の変化点スコアを求める手段を更に備え、
前記演算手段は、
前記第１の変化点スコアと第２の変化点スコアとを加算して変化点スコアを算出するスコア算出手段を更に備える
ことを特徴とする請求項４に記載の変化点予測装置。
時系列的に連続する事象に係る値及び該値に関連する複数の関連事象に係る値を複数の項目毎の項目値として前記値に対応付けて時系列に記録する記録手段を備えるプロセッサが、前記記録手段から前記値又は該値に対応する前記複数の項目毎の項目値を読み出して演算を行ない、記録期間よりも後の予測対象期間における前記事象の変化点を予測する変化点予測方法であって、
前記プロセッサは、
前記記録手段に記録されている複数の項目の項目値と、該項目値に対応する期間から所定期間後の前記事象に係る値とを対応付け、
各項目について対応付けられた前記項目値と前記値との間の相関の強さを示す要因効果値に基づき選択される複数の項目の項目値から、前記事象に係る値の前記予測対象期間における予測値を算出し、
記録期間における前記事象に係る値と、算出された前記予測対象期間における予測値とを連続させた第１の時系列データを作成し、
作成した第１の時系列データに対し、該第１の時系列データのトレンドをトレンドモデルに基づき算出し、
算出したトレンドの勾配を算出し、
算出した勾配のマハラノビスの汎距離を算出し、
算出したマハラノビスの汎距離に基づき前記予測対象期間を含む期間における変化点を検出する
ことを特徴とする変化点予測方法。
時系列的に連続する事象に係る値及び該値に関連する複数の関連事象に係る値を複数の項目毎の項目値として前記値に対応付けて時系列に記録する記録手段を備えるコンピュータに、前記記録手段から前記値又は該値に対応する前記複数の項目毎の項目値を読み出して演算を行ない、記録期間よりも後の予測対象期間における前記事象の変化点を予測させるコンピュータプログラムであって、
前記コンピュータに、
前記記録手段に記録されている複数の項目の項目値と、該項目値に対応する期間から所定期間後の前記事象に係る値とを対応付けるステップと、
各項目について対応付けられた前記項目値と前記値との間の相関の強さを示す要因効果値に基づき選択される複数の項目の項目値から、前記事象に係る値の前記予測対象期間における予測値を算出するステップと、
記録期間における前記事象に係る値と、算出された前記予測対象期間における予測値とを連続させた第１の時系列データを作成するステップと、
作成した第１の時系列データに対し、該第１の時系列データのトレンドをトレンドモデルに基づき算出するステップと、
算出したトレンドの勾配を算出するステップと、
算出した勾配のマハラノビスの汎距離を算出するステップと、
算出したマハラノビスの汎距離に基づき前記予測対象期間を含む期間における変化点を検出するステップと
を実行させることを特徴とするコンピュータプログラム。