(57)【要約】　　　（修正有）【課題】気象庁が発表する長期予報を反映した日々の気象物理量を推定する方法、システム、プログラムおよびプログラムを格納した記憶媒体を提供する。【解決手段】気象物理量を対象とした気象庁発表の長期予報を取得するデータ取得ステップ０２０１、気象物理量の従う時系列モデルを作成する気象時系列モデル作成ステップ０２０２、推定期間における推定地点の気象物理量をシミュレーション回数だけシミュレートする気象物理量シミュレーションステップ０２０３、気象物理量のシミュレーション結果を気象物理量推定結果として出力する気象物理量推定結果出力ステップ０２０４を有する。さらに、気象時系列モデル作成ステップ０２０２を、気象時系列モデルのパラメータ推定ステップ０２０２１および気象時系列モデルのパラメータ補正ステップ０２０２２から構成する。【選択図】　図２
2004117228.tif 000002
【特許請求の範囲】
【請求項１】
過去に観測した気温・降水量などの気象物理量の履歴と、気象物理量を推定する推定期間、推定する推定地点およびシミュレーション回数を含む気象物理量の推定条件と、前記推定地点の前記推定期間における気象物理量を対象とした気象庁発表の長期予報とを入力とし、入力した前記気象物理量の履歴と前記長期予報に基づいて、前記推定地点の前記推定期間における気象物理量に対する気象時系列モデルを作成し、この作成したモデルに基づいて前記シミュレーション回数だけシミュレートした結果を気象物理量推定結果として出力することを特徴とする気象物理量の推定方法。
【請求項２】
過去に観測した気温・降水量などの気象物理量の履歴と、気象物理量を推定する推定期間、推定する推定地点およびシミュレーション回数を含む気象物理量の推定条件と、前記推定地点の前記推定期間における気象物理量を対象とした気象庁発表の長期予報とを取得するデータ取得ステップ、前記データ取得ステップで取得した過去の気象物理量の履歴と長期予報に基づいて、前記推定地点の前記推定期間における気象物理量に対する気象時系列モデルを作成する気象時系列モデル作成ステップ、前記気象時系列モデル作成ステップで作成したモデルに基づいて前記シミュレーション回数だけシミュレートする気象物理量シミュレーションステップ、前記気象物理量シミュレーションステップのシミュレーション結果を気象物理量推定結果として出力する気象物理量推定結果出力ステップ、を備える気象物理量の推定方法であって、
前記気象時系列モデル作成ステップは、前記データ取得ステップで取得した過去の気象物理量の履歴に基づいてモデルのパラメータを推定する気象時系列モデルのパラメータ推定ステップと、前記データ取得ステップで取得した長期予報に基づいて前記気象時系列モデルのパラメータ推定ステップで推定したパラメータに対して補正処理を行う気象時系列モデルのパラメータ補正ステップから成ること、を特徴とする気象物理量の推定方法。
【請求項３】
請求項２に記載の気象物理量の推定方法であって、
前記データ取得ステップで取得する長期予報は、前記推定期間を包含する予報期間と、前記予報期間における期間平均気温に対する「平年並み」の階級の区間および「低い」「平年並み」「高い」の３階級の確率分布を含み、前記気象時系列モデルのパラメータ補正ステップは、モデルのパラメータのうち、期間平均気温の確率分布を特徴付けるパラメータに対して、前記データ取得ステップで取得した長期予報に基づいて補正処理を行う気象時系列モデルのパラメータ補正ステップであること、を特徴とする気象物理量の推定方法。
【請求項４】
請求項２に記載の気象物理量の推定方法であって、
前記データ取得ステップで取得する長期予報は、前記推定期間を包含する予報期間と、前記予報期間における期間累積降水量に対する「平年並み」の階級の区間および「少ない」「平年並み」「多い」の３階級の確率分布を含み、前記気象時系列モデルのパラメータ補正ステップは、モデルのパラメータのうち、期間累積降水量の確率分布を特徴付けるパラメータに対して、前記データ取得ステップで取得した長期予報に基づいて補正処理を行う気象時系列モデルのパラメータ補正ステップであること、を特徴とする気象物理量の推定方法。
【請求項５】
過去に観測した気温・降水量などの気象物理量の履歴と、気象物理量を推定する推定期間、推定する推定地点およびシミュレーション回数を含む気象物理量の推定条件と、前記推定地点の前記推定期間における気象物理量を対象とした気象庁発表の長期予報とを取得するデータ取得手段、前記データ取得手段で取得した過去の気象物理量の履歴に基づいて、推定地点の推定期間における気象物理量に対する気象時系列モデルのパラメータを推定すると共に、前記データ取得ステップで取得した長期予報に基づいて、前記推定された気象時系列モデルのパラメータを補正する気象時系列モデル作成手段、前記気象時系列モデル作成手段で作成したモデルを用いて前記推定期間における前記推定地点の日々の気象物理量を前記シミュレーション回数だけシミュレートする気象物理量シミュレーション手段、及び前記気象物理量シミュレーション手段のシミュレーション結果を気象物理量推定結果として出力する気象物理量推定結果出力手段を備える気象物理量の推定システム。
【請求項６】
請求項５に記載の気象物理量の推定システムにおいて、
前記データ取得手段で取得する長期予報は、前記推定期間を包含する予報期間と、前記予報期間における期間平均気温に対する「平年並み」の階級の区間および「低い」「平年並み」「高い」の３階級の確率分布を含み、前記気象時系列モデル作成手段は、気象時系列モデルのパラメータを補正する際、モデルのパラメータのうち、期間平均気温の確率分布を特徴付けるパラメータに対して、前記データ取得手段で取得した長期予報に基づいて補正処理を行う構成を有することを特徴とする気象物理量の推定システム。
【請求項７】
請求項５に記載の気象物理量の推定システムにおいて、
前記データ取得手段で取得する長期予報は、前記推定期間を包含する予報期間と、前記予報期間における期間累積降水量に対する「平年並み」の階級の区間および「少ない」「平年並み」「多い」の３階級の確率分布を含み、前記気象時系列モデル作成手段は、気象時系列モデルのパラメータを補正する際、モデルのパラメータのうち、期間累積降水量を特徴付けるパラメータに対して、前記データ取得手段で取得した長期予報に基づいて補正処理を行う構成を有することを特徴とする気象物理量の推定システム。
【請求項８】
請求項５に記載の気象物理量の推定システムにおいて、ネットワークを介して情報を送受信する通信手段を備え、前期データ取得手段が前記長期予報と前記推定条件を前記ネットワークを介して取得する構成を有すると共に、前記気象物理量推定結果出力手段が前記気象物理量推定結果を前記ネットワークを介して出力する構成を有することを特徴とする気象物理量の推定システム。
【請求項９】
過去に観測した気温・降水量などの気象物理量の履歴を用いてある地点におけるある期間の気象物理量を推定するため、コンピュータを、
過去に観測した気温・降水量などの気象物理量の履歴と、気象物理量を推定する期間、推定する地点およびシミュレーション回数を含む気象物理量の推定条件と、推定地点の推定期間における気象物理量を対象とした気象庁発表の長期予報とを取得するデータ取得モジュール、
前記データ取得モジュールで取得した過去の気象物理量の履歴に基づいて、推定地点の推定期間における気象物理量に対する気象時系列モデルのパラメータを推定する気象時系列モデルのパラメータ推定モジュール、
前記データ取得ステップで取得した長期予報に基づいて前記気象時系列モデルのパラメータ推定モジュールで推定されたパラメータを補正する気象時系列のパラメータ補正モジュールから成り気象物理量の従う気象時系列モデルを作成する気象時系列モデル作成モジュール、
前記気象時系列モデル作成モジュールで作成したモデルを用いて前記推定期間における前記推定地点の日々の気象物理量を前記シミュレーション回数だけシミュレートする気象物理量シミュレーションモジュール、及び
前記気象物理量シミュレーションモジュールのシミュレーション結果を気象物理量推定結果として出力する気象物理量推定結果出力モジュール、
として機能させることを特徴とする気象物理量の推定プログラム。
【請求項１０】
請求項９に記載の気象物理量の推定プログラムにおいて、
前記データ取得モジュールで取得する長期予報は、前記推定期間を包含する予報期間と、前記予報期間における期間平均気温に対する「平年並み」の階級の区間および「低い」「平年並み」「高い」の３階級の確率分布を含み、前記気象時系列モデルのパラメータ補正モジュールは、モデルのパラメータのうち、期間平均気温の確率分布を特徴付けるパラメータに対して、前記データ取得手段で取得した長期予報に基づいて補正処理を行う構成を有することを特徴とする気象物理量の推定プログラム。
【請求項１１】
請求項９に記載の気象物理量の推定プログラムにおいて、
前記データ取得ステップで取得する長期予報は、前記推定期間を包含する予報期間と、前記予報期間における期間累積降水量に対する「平年並み」の階級の区間および「少ない」「平年並み」「多い」の３階級の確率分布を含み、前記気象時系列モデルのパラメータ補正モジュールは、モデルのパラメータのうち、期間累積降水量の確率分布を特徴付けるパラメータに対して、前記データ取得手段で取得した長期予報に基づいて補正処理を行う構成を有することを特徴とする気象物理量の推定プログラム。
【請求項１２】
請求項９記載の気象物理量の推定プログラムにおいて、前記コンピュータを、さらにネットワークを介して情報を送受信する通信プログラムモジュールとして機能させ、前期データ取得モジュールが前記長期予報と前記推定条件を前記ネットワークを介して取得する構成を有すると共に、前記気象物理量推定結果出力モジュールが前記気象物理量推定結果を前記ネットワークを介して出力する構成を有することを特徴とする気象物理量の推定プログラム。
【請求項１３】
過去に観測した気温・降水量などの気象物理量の履歴と、気象物理量を推定する期間、推定する地点およびシミュレーション回数を含む気象物理量の推定条件と、推定地点の推定期間における気象物理量を対象とした気象庁発表の長期予報とを取得するデータ取得モジュール、前記データ取得モジュールで取得した過去の気象物理量の履歴に基づいて、推定地点の推定期間における気象物理量に対する気象時系列モデルのパラメータを推定する気象時系列モデルのパラメータ推定モジュール、前記データ取得ステップで取得した長期予報に基づいて前記気象時系列モデルのパラメータ推定モジュールで推定されたパラメータを補正する気象時系列のパラメータ補正モジュールから成り気象物理量の従う気象時系列モデルを作成する気象時系列モデル作成モジュール、前記気象時系列モデル作成モジュールで作成したモデルを用いて前記推定期間における前記推定地点の日々の気象物理量を前記シミュレーション回数だけシミュレートする気象物理量シミュレーションモジュール、および、前記気象物理量シミュレーションモジュールのシミュレーション結果を気象物理量推定結果として出力する気象物理量推定結果出力モジュールとしてコンピュータを機能させることを特徴とする気象物理量の推定プログラムを格納した記憶媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、気温・降水量などの気象物理量の推定方法に関し、特に過去に観測した気象物理量の履歴と、気象庁が発表する長期予報から気象時系列モデルを作成し、この作成したモデルに基づいて気象物理量を推定する方法に関する。
【０００２】
【従来の技術】
気温、降水量、積雪量などの変動によって被る売上金減少や費用増加をヘッジするための金融派生商品として、天候デリバティブ（取引）が知られている。天候デリバティブ取引は、天候により売上が大きく左右される顧客（企業、団体、法人、組合、…）と、例えば損害保険会社（銀行、…）の間でなされる。気温が平年に比べて低い場合の売上減少リスクをヘッジする場合には、例えば、次のような概要の契約を顧客と損害保険会社が結び、天候デリバティブ取引を行う。
【０００３】
（１）観測期間：○月○日から△月△日まで
（２）観測指標：平均気温（１日の平均気温が××［℃］以下の日数）
（３）観測地点：○○市
（４）行使値：４日（１日の平均気温が××［℃］以下の日が４日以上あった場合に５日目から補償金受け取り）
（５）補償金：１日あたり○○万円（上限額○○万円）
（６）オプション料：○万円
このような契約内容の天候デリバティブ取引を行う場合、観測期間内の観測指標をどのように推定するかによって補償金の額が異なってくる。このため、観測指標の推定は、顧客にとっても損害保険会社などにとっても極めて重要になる。
【０００４】
従来、観測指標となる気温・降水量などの気象物理量の推定には、まず観測地点における過去の気象物理量の履歴から気象時系列モデルを作成し、そのモデルに基づいて観測期間内の気象物理量の変動を、モンテカルロ・シミュレーションなどで推定するという方法が採られてきた。
【０００５】
例として、日々の気温に対する気象時系列モデルで有名なＤｉｓｃｈｅｌモデルを挙げる。Ｄｉｓｃｈｅｌモデルは（数１）のように表される。
【０００６】
【数１】
2004117228.tif 000003
ここでｎはある任意の日（例えば２０００年１月１日）を１日目とした場合のｎ日目（例えばｎ＝３２であれば２０００年２月１日）を表し、Θｎはｎ日目にあたる日付の平年気温（例えばｎ＝３２であれば２月１日の過去数年間の平均気温）、Ｔｎはｎ日目の気温を表している。ここでパラメータβ、μ、σは定数である。またεｎは平均μ、分散σ２の正規分布に従う確率変数である。さらにβは、あるｎ日目の気温Ｔｎが、対象日の過去数年間の平年気温Θｎと、対象日ｎの前日の気温Ｔｎ−１に、それぞれどの程度の割合で依存しているかを表している。３つのパラメータの値は一般に、観測地点における過去数年の気温の履歴から、最小二乗法を用いて推定される。
一旦パラメータが決定したら、初期値Ｔ０と乱数列ε１、ε２、ε３、…を与えれば、次々に気温時系列Ｔ１、Ｔ２、Ｔ３、…を推定する事ができる。推定期間がｎ＝１からｎ＝Ｎならば、Ｔ１、…、ＴＮで一組の推定気温時系列である。このような推定気象時系列を、乱数列を変える事により数千〜数万通り推定する。
そして、各推定気温時系列において損害保険会社が顧客に補償金として支払う金額を計算し、その確率密度分布を計算する。損害保険会社は最終的に、この確率密度分布を元にオプション料を計算する。以上が従来のオプション料計算プロセスの概略である。
また、気温・降水量などの気象物理量の推定は、このようは天候デリバティブ取引を行わない場合も、企業戦略や各種行事の遂行や活動などにとって、重要な意味を持つ。また、気温・降水量などの気象物理量は誰もが関心を持っている事象である。
【０００７】
【発明が解決しようとする課題】
天候デリバティブ取引に代表される経済活動で用いられる気象モデルは、空間的に一点の、すなわち観測地点の気象変動のみに注目した簡潔な時系列モデルである。従来通りモデルのパラメータを過去の気象物理量の履歴のみから推定した場合、そのモデルによる将来予測は観測地点における平均的な傾向に過ぎない。つまり、今年の夏は冷夏になる可能性が高い、といった具体的な年毎の傾向を予想することはできない。
そもそも、そのような詳細な気象現象を予想するためには、気象庁が行っているような、極めて詳細な気象力学モデルに基づいた地球的スケールのシミュレーションを実行するなど、莫大な工数および各種データを解釈する卓越した技能が要求される。
実際気象庁は、そのような大規模なシミュレーション結果などの資料に基づいて、気温と降水量について向こう１か月、３か月の長期予報を発表している。その発表内容は、予報期間の期間平均気温（期間降水量）が「低い（少ない）」「平年並み」「高い（多い）」の３つの階級で予報される。各階級の幅は、１９７１年〜２０００年の３０年間における各階級の出現率が等分（それぞれ３３％）となるように決められている。例えば、２００２年８月１６日に発表された「関東甲信地方　１か月予報（８月１７日から９月１６日までの天候見通し）」では、向こう１か月間の平均気温は「低い確率２０％、平年並みの確率５０％、高い確率３０％」であった。長期予報の精度は年々向上しており、今後もコンピュータの能力向上とともに改善されて行くと期待されている。
【０００８】
しかしながら、長期予報は期間平均気温・期間降水量の確率分布というマクロな情報を提供するだけであり、これに対し天候デリバティブ取引に代表される経済活動では、日々の気温・降水量などの確率分布というミクロな情報が必要とされる場合も多い。実際、Ｄｉｓｃｈｅｌモデルなどの気象時系列モデルは、年毎の傾向を予測することができないという短所があるものの、日々の気温・降水量などの確率分布というミクロな情報を取り扱うことができるという長所があったため、天候デリバティブ取引において普及したと考えられる。
【０００９】
発明者は、もしマクロな長期予報を、ミクロな気象時系列モデルに反映させることができれば、両者の長所を兼ね備えた理想的な気象時系列モデルを作成することができることに注目した。
本特許の課題は、日々の気温・降水量などの気象物理量が経済活動に与える影響をより正確に評価する為に、気象庁が発表する長期予報を反映した日々の気象物理量を推定する方法、システム、プログラムおよびプログラムを格納した記憶媒体を提供することである。
【００１０】
【課題を解決するための手段】
前記課題を解決するために、本発明では、過去に観測した気温・降水量などの気象物理量の履歴、気象物理量を推定する期間、推定する地点およびシミュレーション回数を含む気象物理量推定条件、推定地点の推定期間における気象物理量を対象とした気象庁発表の長期予報を取得するデータ取得ステップ、前記データ取得ステップで取得した過去の気象物理量の履歴と長期予報に基づいて、推定地点の推定期間における気象物理量に対する気象時系列モデルを作成しする気象時系列モデル作成ステップ、前記気象時系列モデル作成ステップで作成したモデルに基づいてシミュレーション回数だけシミュレートする気象物理量シミュレーションステップ、前記気象物理量シミュレーションステップのシミュレーション結果を気象物理量推定結果として出力する気象物理量推定結果出力ステップから気象物理量の推定方法を構成する。さらに、前記気象時系列モデル作成ステップを、前記データ取得ステップで取得した過去の気象物理量の履歴に基づいてモデルのパラメータを推定する気象時系列モデルのパラメータ推定ステップと、前記データ取得ステップで取得した長期予報に基づいて前記気象時系列モデルのパラメータ推定ステップで推定したパラメータに対して補正処理を行う気象時系列モデルのパラメータ補正ステップから構成する。
【００１１】
【発明の実施の形態】
以下、図面を参照して本発明の実施形態を詳細に説明する。
実施例として、先ほど述べたＤｉｓｃｈｅｌモデルを取り上げる。まずその理論的側面を説明し、次に理論に基づいた実践的側面について説明する。
＜理論的側面＞
はじめにＤｉｓｃｈｅｌモデルによる期間平均気温の分布は正規分布になることを示し、次に長期予報の確率分布を正規分布に変換する方法を示す。さらに前者の正規分布を、後者の正規分布に一致させるようにパラメータに補正処理を行う方法と、その物理的原理について述べる。
【００１２】
まず、Ｄｉｓｃｈｅｌモデルの生成する気温時系列の期間平均気温の分布が正規分布になることを示す。まず、（数１）のｎ日目の気温Ｔｎは、初期値Ｔ０と平年気温列Θ１、Θ２、…を用いて次のように下すことができる。
【００１３】
【数２】
2004117228.tif 000004
第一項目と第二項目は確率的な変動には依存しない項であり、確率的な変動は第三項に集約されている。さらに、長期予報の期間平均気温に対応する物理量として、１日目からＮ日目までのＮ日間の期間平均気温τを次のように定義する。
【００１４】
【数３】
2004117228.tif 000005
τは次のように書き下すことができる。
【００１５】
【数４】
2004117228.tif 000006
ここでも第一項目と第二項目は確率的な変動には依存しない項であり、確率的な変動は第三項に集約されている。εｉは平均μ分散σ２の正規分布に従っているので、結局τは平均Ｅ［τ］、分散Ｖ［τ］の正規分布に従っていることがわかる。ただしＥ［τ］とＶ［τ］は次式で与えられる。
【００１６】
【数５】
2004117228.tif 000007
以上のようにＤｉｓｃｈｅｌモデルの生成する気温時系列の期間平均気温の分布が正規分布になることが示された。
一方、長期予報はτが従う確率分布を「低い確率」、「平年並みの確率」、「高い確率」で予報しており、（数５）で特徴付けられる正規分布と比較するのに適していない。
そこで、長期予報の確率分布を正規分布に変換する方法について述べる。「低い確率」をｐ１、「高い確率」をｐ２と書くことにする。この時「平年並みの確率」は当然１−ｐ１−ｐ２である。ここで、平均０分散１の標準正規分布に従っている確率変数ｘを考え、ｘ１、ｘ２を次のように定義する。
【００１７】
【数６】
2004117228.tif 000008
与えられたｐ１、ｐ２に対して、（数６）の方程式をそれぞれｘ１、ｘ２について数値的に解くのは容易である。実際、多くの数値計算ソフトにそのような機能が付属しているので、その解法については言及しない。ここで、長期予報によって与えられた平年並みの階級の区間を［τ１，τ２］と書くことにする。この区間は実際には、例えば「向こう１か月間の平均気温は、平均２５．８℃、気温平年差−０．３℃〜＋０．４℃」という形式で発表される。つまりこの場合、平年並みの階級の区間は［２５．５℃，　２６．２℃］である。さらに、次のような確率変数ｙを考える。
【００１８】
【数７】
2004117228.tif 000009
ただしＣ１、Ｃ２はこれから決定する定数である。明らかにｙは平均Ｃ１分散Ｃ２２の正規分布に従う。ｙ１＝Ｃ２ｘ１＋Ｃ１、ｙ２＝Ｃ２ｘ２＋Ｃ１と書くことにしよう。この時、ｙ１＝τ１、ｙ２＝τ２となるようにＣ１、Ｃ２を決定すれば、ｙは次のように各階級の確率分布が長期予報と等しい正規分布に従う。
【００１９】
【数８】
2004117228.tif 000010
具体的に連立方程式ｙ１＝Ｃ２ｘ１＋Ｃ１、ｙ２＝Ｃ２ｘ２＋Ｃ１をＣ１、Ｃ２について解くと次式を得る。
【００２０】
【数９】
2004117228.tif 000011
このようにして長期予報を、各階級の確率分布が長期予報のそれと等しい正規分布に変換することができる。以後、この変換された分布を、長期予報の正規分布と呼ぶ事にする。
【００２１】
二つの正規分布を比較するには、平均と分散を比較すれば良い。Ｄｉｓｃｈｅｌモデルの正規分布の平均と分散はＥ［τ］とＶ［τ］であり、（数５）で与えられる。また長期予報の正規分布の平均と分散はＣ１、Ｃ２２であり、（数９）で与えられる。両者は当然一般には一致しない。
しかしながら、使用されている情報量と気象モデルの精巧さなどを考慮すれば、長期予報の方が未来の気象現象をより的確に予測していると考えられる。
そこで、Ｄｉｓｃｈｅｌモデルの正規分布が長期予報の正規分布に一致するように、パラメータを決定する方法について以下に述べる。
長期予報だけから全てのパラメータが決定できれば一番良いのだが、Ｄｉｓｃｈｅｌモデルの場合には、それだけでは情報が足りない。（数５）、（数９）から平均と分散がそれぞれ等しいという連立方程式を立てると（数１０）、（数１１）のようになる。
【００２２】
【数１０】
2004117228.tif 000012
【数１１】
2004117228.tif 000013
方程式２本に対して求めるべきパラメータ（β、μ、σ）が３つあるので、解を一意的に決定できないことがわかる。つまり、新たな情報が必要である。
【００２３】
そこで、過去の気温履歴から有益な情報を引き出そう。例として観測地点東京における１９７１年１月１日から２０００年１２月３１日まで日毎の平均気温データに注目する。このデータを３か月毎（季節毎）の区間に分け、区間毎にパラメータを推定した結果を図１に示す。各パラメータは区間毎に大きく揺らいでいることがわかる。このとき各パラメータの平均と標準偏差は表１のような結果であった。
【００２４】
【表１】
2004117228.tif 000014
この揺らぎの主な原因の一つは、例えば「今年の夏は猛暑である」あるいは「今年の冬は暖冬である」といった平年よりも偏った気象現象の発現である。逆に、パラメータが表１の標準偏差程度（Δβ、Δμ、Δσと書く）ゆらいだ場合、（数５）で与えられるＥ［τ］とＶ［τ］はどれだけ変化するだろうか。具体的に評価すると（数１２）のような結果を得る。
【００２５】
【数１２】
2004117228.tif 000015
但し、Ｅ［τ］とＶ［τ］は２００１年６月１日から２００１年８月３１日までの期間平均気温τについて計算した値である。また、各パラメータは１９７１年１月１日から２０００年１２月３１日までのデータから推定し、Ｔ０は２００１年５月３１日の一日平均気温とした。この結果は次のように解釈することができる。仮にτの平均値（すなわちＥ）が１℃変化したとしよう。Ｅの変化はβのゆらぎが原因であるとすると、βのゆらぎ幅はΔβの４倍程度必要になるが、標準偏差の４倍程度の大きなゆらぎが起きる確率は１万分の１より小さい。これに対し、Ｅの変化はμのゆらぎが原因であるとすると、μのゆらぎ幅はほぼΔμでよく、そのようなゆらぎが起きる確率は約３０％である。従って、Ｅの変化を説明するのに最も確実性が高い。このようにしてＥの変化を回帰する先はμであると結論できる。同様にして、Ｖの変化を回帰する先はβであると結論できる。還元すれば、σのゆらぎは近似的に無視でき、過去の気温履歴から推定した値を利用しても差し支えない。
この物理的性質に着目し、まず一旦過去の気温履歴からβ、μ、σを最小二乗法によって決定する。次に長期予報に基づいて、（数１１）を満たすようにβを数値的に解き直し、その補正処理後のβを（数１０）に代入して、μについて解き直す。このようにして、補正処理後のパラメータセットを得ることができる。
以上のようにして、長期予報を気象時系列モデルのパラメータに反映させて、今年の夏は冷夏になる可能性が高い、といった空間的に広大な領域の物理情報を実効的に取り込むことが可能となった。ただし、ここで展開した議論はＤｉｓｃｈｅｌモデルに限定されたものではなく、幾つか提案されている他の気象時系列モデルであっても同様の議論が可能である。
＜実践的側面＞
図２は本発明に係わる気象物理量推定方法の実施形態を表すフローチャートである。
図２に示すように、本実施形態の気象物理量推定方法は、過去に観測した気象物理量の履歴情報である気象物理量履歴０２１１および、気象庁発表の長期予報０２１２および、気象物理量を推定する条件である推定条件０２１３を取得するデータ取得ステップ０２０１、データ取得ステップ０２０１にて取得した情報をもちいて気象物理量の時系列モデルを作成する気象時系列モデル作成ステップ０２０２、気象時系列モデル作成ステップ０２０２にて作成した気象系列モデルを用いて、気象物理量をシミュレートする気象物理量シミュレーションステップ０２０３、気象物理量シミュレーションステップ０２０３のシミュレーション結果を気象物理量推定結果０２１４として出力する気象物理量推定結果出力ステップ０２０４から成る。さらに、気象時系列モデル作成ステップ０２０２は、気象物理量履歴０２１１から気象時系列モデルのパラメータの値を推定する気象時系列モデルのパラメータ推定ステップ０２０２１および、気象時系列モデルのパラメータ推定ステップで推定したパラメータに対し長期予報０２１２から補正処理を行う気象時系列モデルのパラメータ補正ステップ０２０２２から成る。
次に、図２中の気象物理量履歴０２１１の一例を、図３を用いて説明する。
本実施形態の気象物理量履歴０２１１は、日付０３１および各観測地点における各日の１日平均気温０３２からなる。図３において、例えば１９７２年０１月０１日における東京の１日平均気温は４．１［℃］である。
次に、図２中の長期予報０２１２の一例を、図４を用いて説明する。
本実施形態の長期予報０２１２は、予報期間０４１、平年並みの範囲０４２、および「低い」「平年並み」「高い」の３階級の確率分布０４３から成る。図４は、気象庁が発表する長期予報の発表内容を表している。
次に、図２中の推定条件０２１３の一例を、図５を用いて説明する。
本実施形態の推定条件０２１３は、１日平均気温の推定期間０５１、推定地点０５２およびシミュレーション回数０５３から成る。図５は、東京における０６月０１日から０８月３１日までの１日平均気温時系列をシミュレーションによって１００００通り生成することを表している。
次に、図２中のデータ取得ステップ０２０１を、図３、５、１０を用いて説明する。
図５に示す推定条件０２１３では推定地点０５２が東京のため、図２中のデータ取得ステップ０２０１では、図３に示す気象物理量履歴０２１１から取得する気象物理量履歴は、図１０に示すように、東京における過去の気象物理量となる。ただし、本実施形態のおいては、図４中の予報期間０４１と、図５中の推定期間０５１とは、予報期間０４１が推定期間０５１と一致する関係、あるいは予報期間０４１が推定期間０５１を含む関係にあるものとする。
次に本実施形態の気象時系列モデルのパラメータ推定ステップ０２０２１の一例を、図６を用いて説明する。この気象時系列モデルのパラメータ推定ステップ０２０２１は、過去の気象データから気象時系列モデルのパラメータを推定するものである。
本実施形態の気象時系列モデルのパラメータ推定ステップ０２０２１では、（数１）で与えられるＤｉｓｃｈｅｌモデルに関するパラメータβ、μ、σを推定する。
図６に示すように、図２の気象時系列モデルのパラメータ推定ステップ０２０２１は、初期化ステップ０６０１、観測気温配列作成ステップ０６０２、平年気温配列作成ステップ０６０３、観測点配列作成ステップ０６０４、回帰分析ステップ０６０５から成る。
初期化ステップ０６０１では、平均気温履歴０２１１（図１０）を参照し、履歴の総年数Ｙを格納する。例えば図１０は３１年分のデータなのでＹ＝３１である。観測気温配列作成ステップ０６０２では、平均気温に履歴０２１１（図１０）を参照し、うるう年の場合には２月２９日をスキップしながら、観測気温配列　Ｔ［Ｉ］（Ｉ＝１、…、３６５＊Ｙ）を作成する。平年気温配列作成ステップでは、観測気温配列からＹ年間の同じ日付の日についての気温の平均、すなわち平年気温配列　Θ［Ｉ］（Ｉ＝１、…、３６５）を、（数１３）によって作成する。
【００２６】
【数１３】
2004117228.tif 000016
ここで、（数１）に対して（数１４）のような回帰式を立てる。
【００２７】
【数１４】
2004117228.tif 000017
ただしδは平均０の残差である。この回帰式に対して、説明変数ｘと目的変数ｙをそれぞれ、（数１５）としたＭ＝３６５＊Ｙ−１個の観測点を用意する。
【００２８】
【数１５】
2004117228.tif 000018
そして、これらの観測点によって回帰分析を行い、βとμおよび残差δの分散であるσ２　　　　　を求める。その為に、観測点配列作成ステップでは、観測気温配列Ｔ［Ｉ］と平年気温配列Θ［Ｉ］から、目的変数配列Ｙ［Ｉ］（Ｉ＝１、…、Ｍ）と説明変数配列Ｘ［Ｉ］（Ｉ＝１、…、Ｍ）を（数１６）に従って作成する。
【００２９】
【数１６】
2004117228.tif 000019
回帰分析ステップでは、観測点配列作成ステップで作成したＸ［Ｉ］とＹ［Ｉ］をもとに、最小二乗法によって導出した（数１７）に従ってβ、μ、σ２を推定する。
【００３０】
【数１７】
2004117228.tif 000020
次に本実施形態の気象時系列モデルのパラメータ補正ステップ０２０２２の一例を、図７を用いて説明する。この気象時系列モデルのパラメータ補正ステップ０２０２１は、過去の気象データから決定したパラメータに対し、長期予報の情報を取り込むように補正処理を行うものである。
図７に示すように、図２の気象時系列モデルのパラメータ補正ステップ０２０２２は、長期予報データ格納ステップ０７０１、推定条件格納ステップ０７０２、ｘ１およびｘ２の計算ステップ０７０３、β計算ステップ０７０４、μ計算ステップ０７０５から成る。
【００３１】
長期予報データ格納ステップ０７０１では、図２のデータ取得ステップで取得した長期予報から、「低い確率」と「高い確率」をそれぞれｐ１とｐ２に、平年並みの区間の下限と上限をそれぞれτ１とτ２に格納する。また、予報期間の日数を計算してＮに格納し、平年気温配列Θ［Ｉ］（Ｉ＝１、…、３６５）から予報期間の平年気温配列Θ’［Ｉ］（Ｉ＝１、…、Ｎ）を格納し、平年気温配列Θ［Ｉ］（Ｉ＝１、…、３６５）から予報期間の前日の気温をＴ０として格納する。正規分布の逆関数計算ステップ０７０２では、方程式（数６）をｘ１とｘ２について数値的に解く。β計算ステップ０７０３では、方程式（数１１）をβについて数値的に解く。ただし、σ２は図６の回帰分析ステップ０６０５で取得した値を用いる。また、方程式（数１１）をβについて数値的に解くのは、ニュートン・ラプソン法などを用いれば容易に解くことができるので詳細な解法については言及しない。μ計算ステップ０７０４では、（数１２）の方程式をμについて解いた（数１８）によってμを計算する。ただし、βはβ計算ステップ０７０３で取得した値を用いる。
【００３２】
【数１８】
2004117228.tif 000021
次に本実施形態の気象物理量シミュレーションステップ０２０３の一例を、図８を用いて説明する。この気象物理量シミュレーションステップ０２０３は、気象時系列モデル作成ステップ０２０２で作成した気象時系列モデルに基づいてシミュレーションを実行するものである。
【００３３】
図８に示すように、図２の気象物理量シミュレーションステップ０２０３は、推定条件格納ステップ０８０１、シミュレーション初期化ステップ０８０２、乱数発生ステップ０８０３、気温推定ステップ０８０４、推定期間終了判定ステップ０８０５、推定日数加算ステップ０８０６、推定気温時系列蓄積ステップ０８０７、シミュレーション回数終了判定ステップ０８０８、シミュレーション番号加算ステップ０８０９から成る。
【００３４】
推定条件格納ステップ０８０１では、推定期間の日数を計算してＮ’に格納し、シミュレーション回数をＬに格納し、シミュレーション番号Ｐを１に初期化する。シミュレーション初期化ステップ０８０２では、推定気温時系列Ｓ［Ｉ］（Ｉ＝０、…、Ｎ’）を、Ｓ［０］＝Ｔ０、Ｓ［Ｉ］＝０（Ｉ＝１、…、Ｎ’）と初期化し、推定日数カウンタＣを１に初期化する。乱数発生ステップ０８０３では、平均μ分散σ２の正規分布に従った乱数Ｕを発生する。ただし、μは図７のμ計算ステップ０７０５で取得した値を、またσ２は図６の回帰分析ステップ０６０５で取得した値を用いる。気温推定ステップ０８０４では、（数１９）に従ってＣ日目の気温Ｓ［Ｃ］を推定する。
【００３５】
【数１９】
2004117228.tif 000022
推定期間終了判定ステップ０８０５では、推定日数カウンタＣが推定期間日数Ｎ’に等しいかどうかを判定し、等しくなければ推定日数加算ステップ０８０６に、また等しければ推定気温時系列蓄積ステップ０８０７に、それぞれ処理が移る。推定日数加算ステップ０８０６では、推定日数Ｃをカウントアップする。推定気温時系列蓄積ステップ０８０７では、推定気温時系列Ｓ［Ｉ］（Ｉ＝０、…、Ｎ’）をデータ蓄積用二次元配列Ｒ［Ｊ］［Ｉ］（Ｊ＝１、…、Ｌ；Ｉ＝１、…、Ｎ’）に、Ｒ［Ｐ］［Ｉ］＝Ｓ［Ｉ］（Ｉ＝１、…、Ｎ’）となるようにコピーする。シミュレーション終了判定ステップ０８０８ではシミュレーション番号Ｐがシミュレーション回数Ｌに等しいかどうかを判定し、等しくなければシミュレーション番号加算ステップ０８０９に、また等しければ終了し気象物理量推定結果出力ステップ０２０４に、それぞれ処理が移る．
次に、図２の気象物理量推定結果出力ステップ０２０４を、図２、８を用いて説明する。図２の気象物理量推定結果出力ステップ０２０４は、図８の推定気温時系列蓄積ステップ０８０７において蓄積されたデータ蓄積用二次元配列Ｒ［Ｊ］［Ｉ］（Ｊ＝１、…、Ｌ；Ｉ＝１、…、Ｎ’）などを元に、図２の気象物理量推定結果０２１４を出力する。
次に、図２の気象物理量推定結果０２１４の一例を、図９を用いて説明する。本実施形態の気象物理量推定結果０２１４は、日付０９１、シミュレーション番号０９２、および各シミュレーション番号の各日付における推定気温０９３からなる。図９において、例えばシミュレーション番号２の０７月０１日の１日推定気温は２５．３［℃］である。
次に、本発明に係わる気象物理量推定システムの実施形態（第一実施形態）を、図１１から図１４を用いて説明する。
本実施形態の気象物理量推定システムは、図１１のように表示装置１１１、ＣＰＵ（演算処理装置）１１２、システムの起動プログラム１１３１を格納したＲＯＭ（Ｒｅａｄ　Ｏｎｌｙ　Ｍｅｍｏｒｙ）１１３、キーボード装置１１４、マウス装置１１５、主記憶装置１１６、オペレーティングシステム１１７１、気象物理量推定プログラム１１７２および気象物理量履歴０２１１を格納した外部記憶装置１１７が通信バス１１８で繋がる構成を取る。
【００３６】
まず物理量推定システムのハードウェアを起動すると、図１２のようにオペレーティングシステム１１７１を主記憶装置１１６へロードする。さらに、気象物理量推定システムのアプリケーションを起動すると、図１３のようにオペレーティングシステム１１７１および気象物理量推定プログラム１１７２を主記憶装置１１６へロードする。さらに、気象物理量の推定を実行すると、図１４のように気象物理量推定プログラム１１７１が気象物理量履歴０２１１、並びに表示装置１１１、キーボード装置１１４、マウス装置１１５を通じて入力された長期予報０２１２と推定条件０２１３を読み込んで気象物理量の推定を実行し、気象物理量推定結果０２１４を外部記憶装置１１７へ出力する。
【００３７】
次に図１１から図１４の気象物理量推定プログラム１１７２の一例を、図１５を用いて説明する。
【００３８】
本実施形態の気象物理量推定プログラムは、データ取得モジュール１５１、気象時系列モデル作成モジュール１５２、気象物理量シミュレーションモジュール１５３および気象物理量推定結果出力モジュール１５４から成る。さらに、気象時系列モデル作成モジュール１５２は気象時系列モデルのパラメータ推定モジュール１５２１および気象時系列モデルのパラメータ補正モジュール１５２２から成る。
データ取得モジュール１５１、気象時系列モデル作成モジュール１５２、気象物理量シミュレーションモジュール１５３、気象物理量推定結果出力モジュール１５４、気象時系列モデルのパラメータ推定モジュール１５２１、気象時系列モデルのパラメータ補正モジュール１５２２は、それぞれ図２中のデータ取得ステップ０２０１、気象時系列モデル作成ステップ０２０２、気象物理量シミュレーションステップ０２０３、気象物理量推定結果出力ステップ０２０４、気象時系列モデルのパラメータ推定ステップ０２０２１、気象時系列モデルのパラメータ補正ステップ０２０２２に従った処理を行うプログラムモジュールである。
【００３９】
第二実施形態の気象物理量推定システムの実施形態（第二実施形態）を、図１６から図１９を用いて説明する。なお、第二実施形態は、通信ネットワークを介して情報の送受信を行う実施形態である。
【００４０】
第二実施形態の気象物理量推定システムは、図１６のようにクライアント１６１、通信ネットワーク１６２、サーバ１６３、クライアント１６１と通信ネットワーク１６２を接続する通信線１６４およびサーバ１６３と通信ネットワーク１６２を接続する通信線１６５から成る。さらに、クライアント１６１は表示装置１６１１、ＣＰＵ１６１２、システムの起動プログラム１６１３１を格納したＲＯＭ１６１３、キーボード装置１６１４、マウス装置１６１５、主記憶装置１６１６、オペレーティングシステム１６１７１、通信モジュール１６１７２を格納した外部記憶装置１６１７、通信装置１６１８が通信バス１６１９で繋がる構成を取る。ここで通信モジュール１６１７２は通信装置１６１８および通信線１６４を介してクライアント１６１と通信ネットワーク１６２の情報をやり取りするプログラムモジュールである。また、サーバ１６３はＣＰＵ１６３１、システムの起動プログラム１６３２１を格納したＲＯＭ１６３２、通信装置１６３３、主記憶装置１６３４、オペレーティングシステム１６３５１、気象物理量推定プログラム１６３５２、通信モジュール１６３５３および気象物理量履歴０２１１を格納した外部記憶装置１６３５が通信バス１６３６で繋がる構成を取る。ここで通信モジュール１６３５３は通信装置１６３３および通信線１６５を介してサーバ１６３と通信ネットワーク１６２の情報をやり取りするプログラムモジュールである。
【００４１】
気象物理量推定システムのクライアントおよびサーバのハードウェアを起動すると、図１７のようにクライアント１６１側ではオペレーティングシステム１６１７１を主記憶装置１６１６へロードし、サーバ１６３側ではオペレーティングシステム１６３５１を主記憶装置１６３４へロードする。さらに、気象物理量推定システムのクライアントおよびサーバのアプリケーションを起動すると、図１８のようにクライアント１６１側では通信モジュール１６１７２を主記憶装置１６１６へロードし、サーバ１６３側では気象物理量推定プログラム１６３５２及び通信モジュール１６３５３を主記憶装置１６３４へロードする。
【００４２】
さらに、気象物理量の推定を実行すると、図１９のように、気象物理量推定プログラム１６３５２が気象物理量履歴０２１１、並びにクライアントの表示装置１６１１、キーボード装置１６１４、マウス装置１６１５を通じて入力された長期予報０２１２と推定条件０２１３を読み込んで気象物理量の推定を実行し、気象物理量推定結果０２１４をクライアントの外部記憶装置１６１７へ出力する。なお、気象物理量推定プログラム１６３５２は図１５で示した実施形態と同様の構成をとる。
次に本発明に係わる第一・第二実施形態の気象物理量推定システムにおける長期予報の入力画面の一例を、図２０を用いて説明する。
【００４３】
一例として示す長期予報の入力画面は、予報期間入力プルダウンメニュー２０１、確率分布入力プルダウンメニュー２０２から成る。本入力画面にて入力された結果が長期予報０２１２として気象物理量推定プログラムに入力される。なお、長期予報の入力画面に表示する情報はこの限りではない。ちなみに、第二実施形態では、入力画面表示用のデータは、例えばサーバ１６３から送信され、クライアント１６１の図示していないブラウザプログラムにより入力画面として構築される。
次に本発明に係わる第一・第二実施形態の気象物理量推定システムにおける推定条件の入力画面の一例を、図２１を用いて説明する。
【００４４】
一例として示す推定条件の入力画面は、推定期間入力プルダウンメニュー２１１、推定地点入力２１２およびシミュレーション回数入力プルダウンメニュー２１３から成る。本入力画面にて入力された結果が推定条件０２１３として気象物理量推定プログラムに入力される。なお、推定条件の入力画面に表示する情報はこの限りではない。ちなみに、第二実施形態では、入力画面表示用のデータは、例えばサーバ１６３から送信され、クライアント１６１の図示していないブラウザプログラムにより入力画面として構築される。
以上説明した本発明は、前記した実施形態に限定されることなく幅広く変形実施することができる。
まず、ここで挙げたＤｉｓｃｈｅｌモデルは一例に過ぎず、他に提案されている気象時系列モデルであっても同様な方法でモデルに長期予報を取り込むことができる。また、例えば１か月長期予報では、向こう１か月間の期間平均気温・期間累積降水量の確率分布の他に、向こう１か月の第１週・第２週・第３〜４週毎の期間平均気温・期間累積降水量の確率分布を提供しているが、これらの情報を取り込むように本実施形態を拡張することは容易である。例えば、第１週・第２週・第３〜４週毎にパラメータ補正処理を行い、シミュレーションを行う際各区間に対応したパラメータを用いればよい。
また、例えば、通信ネットワークを介した第二実施形態は、クライアントに対して課金する有償サービスにすることもできる。また、課金しない有償サービスにすることもできる。
【００４５】
【発明の効果】
本発明によれば、将来の日々の気温・降水量などの気象物理量を、過去の履歴にのみ基づくのでなく、長期予報を反映して推定できるので、経済活動に与える影響をより正確に評価することができる。
【図面の簡単な説明】
【図１】パラメータ推定値の３か月（季節毎）の推移を示すグラフである。
【図２】本発明に係わる実施形態の気象物理量の推定方法を示すフローチャートである。
【図３】図２の気象物理量履歴の一例を示すテーブルである。
【図４】図２の長期予報の一例を示すテーブルである。
【図５】図２の推定条件の一例を示すテーブルである。
【図６】図２の気象時系列モデルのパラメータ推定ステップの一例を示すフローチャートである。
【図７】図２の気象時系列モデルのパラメータ補正ステップの一例を示すフローチャートである。
【図８】図２の気象物理量シミュレーションステップの一例を示すフローチャートである。
【図９】図２の気象物理量推定結果の一例を示すテーブルである。
【図１０】図５の推定条件に基づいて図３から取得した気象物理量履歴の一例を示したテーブルである。
【図１１】本発明に係わる第一実施形態の気象物理量推定システムの構成例を示すブロック図である。
【図１２】図１１の気象物理量推定システムにおけるハードウェア起動時の状態の例を示すブロック図である。
【図１３】図１１の気象物理量推定システムにおけるアプリケーション起動時の状態の例を示すブロック図である。
【図１４】図１１の気象物理量推定システムにおける気象物理量推定実行時の状態の例を示すブロック図である。
【図１５】図１１の気象物理量推定システムにおける気象物理量推定プログラムモジュール一例を示すブロック図である。
【図１６】本発明に係わる第二実施形態の気象物理量推定システムの構成例を示すブロック図である。
【図１７】図１６の気象物理量推定システムにおけるハードウェア起動時の状態の例を示すブロック図である。
【図１８】図１６の気象物理量推定システムにおけるアプリケーション起動時の状態の例を示すブロック図である。
【図１９】図１６の気象物理量推定システムにおける気象物理量推定実行時の状態の例を示すブロック図である。
【図２０】図１１・図１６の気象物理量推定システムにおける長期予報入力画面の一例を示す図である。
【図２１】図１１・図１６の気象物理量推定システムにおける推定条件入力画面の一例を示す図である。
【符号の説明】
０２１１…気象物理量履歴、０２１２…長期予報、０２１３…推定条件、０２０１…データ取得ステップ、０２０２…気象時系列モデル作成ステップ、０２０３気象物理量シミュレーションステップ、０２０４…気象物理量推定結果出力ステップ、０２１４…気象物理量推定結果、０２０２１…気象時系列モデルのパラメータ推定ステップ、０２０２２…気象時系列モデルのパラメータ補正ステップ。
【図１】
2004117228.tif 000023
【図２】
2004117228.tif 000024
【図３】
2004117228.tif 000025
【図４】
2004117228.tif 000026
【図５】
2004117228.tif 000027
【図６】
2004117228.tif 000028
【図７】
2004117228.tif 000029
【図８】
2004117228.tif 000030
【図９】
2004117228.tif 000031
【図１０】
2004117228.tif 000032
【図１１】
2004117228.tif 000033
【図１２】
2004117228.tif 000034
【図１３】
2004117228.tif 000035
【図１４】
2004117228.tif 000036
【図１５】
2004117228.tif 000037
【図１６】
2004117228.tif 000038
【図１７】
2004117228.tif 000039
【図１８】
2004117228.tif 000040
【図１９】
2004117228.tif 000041
【図２０】
2004117228.tif 000042
【図２１】
2004117228.tif 000043