第8章：エネルギー貯蔵 I

章の概要

本章では、非常に変動の激しい価格を持つ送電網から電力を購入したり売却したりする際に発生する、一見かなり単純に見える在庫問題を対象とする、エネルギー貯蔵の問題を考察する。最初の6つの章とは対照的に、本章ではこうした確率的価格を記述するために、より豊富なモデル群を用いる。これにより、不確実性をモデル化する際に直面しうる複雑さの一端が垣間見えることになる。本章では、古典的な時系列モデル、(価格スパイクを捉えるための)ジャンプ拡散モデル、分位点分布、さらには分位点分布と標準正規分布を組み合わせたハイブリッドモデルなど、価格過程のための様々なモデルを一通り紹介し、正規性に依存する手法を利用する道を開く。

続いて、一連の方策について説明する。まず基本的な「安く買って高く売る」方策(方策関数近似の一形態)から始め、次に第5章で最初に見たベルマン方程式を近似する手法へと移行する。まず、ほとんどすべての問題にとって計算上実行不可能な、ベルマン方程式の素朴な記述から出発し、その後、後方近似動的計画法(ADP)、前方ADP、そしてパラメータ調整と組み合わせた前方ADPを用いるハイブリッド戦略として知られる各種のバリエーションを一通り紹介する。

ナラティブ

ニュージャージー州は、3,500メガワット(MW)の洋上風力発電を開発しようとしている。課題の一つは、風力(特に洋上風力)が非常に変動しやすいことである。この変動が送電網に与える影響は、(中間的な範囲において)風力発電の出力が風速の3乗に比例して増加するという性質によって増幅される。この変動性は図8.1に示されている。

風力によるエネルギーは、米国中西部、欧州沖の沿岸地域、ブラジル北東部、中国北部地域(ほんの一例を挙げただけだが)など、風の強い地域で普及してきた。地域社会(や企業)は、炭素排出量の削減や送電網への依存の最小化を目的として、再生可能エネルギー(風力や太陽光)に投資することがある。

しかし、こうしたプロジェクトによって地域社会が送電網をポートフォリオから完全に排除できることは、極めて稀である。一般的な慣行としては、再生可能エネルギー源(風力や太陽光)は送電網に直接売却させる一方、企業は送電網から購入するという形が取られる。これはヘッジとして有用な場合がある。というのも、企業は価格スパイクの際(価格は1メガワット時(mwh)あたり20ドルから300ドル以上に跳ね上がることがある)に多額の利益を得られ、これがそうした期間における電力購入コストを相殺するからである。

再生可能エネルギーにおける主要な難点は、変動性への対処である。一つの解決策は、再生可能エネルギー源からのエネルギーを単純に送電網に流し込み、送電網の容量を利用してこの変動性を吸収するというものだが、変動の山と谷を平滑化するために貯蔵(特に電池による貯蔵)を用いることへの関心も高まっている。再生可能エネルギー源の変動性を平滑化することに加え、価格スパイクを利用するために電池を用いること、すなわち電力が安いとき(価格はマイナスになることさえある)に購入し、高いときに売り戻すことへの関心も高まっている。送電網上の電力価格の変動性を利用して、安いときに買い、高いときに売ることは、電池アービトラージとして知られている。

本章では、図8.2に示す構成を用いて、エネルギー貯蔵における多くのモデリングおよびアルゴリズム上の課題を説明する。この問題は、以下のようなあらゆる在庫/貯蔵問題に対する洞察を与えてくれる。

投資信託における現金保有 – 銀行は、投資資本のうちどれだけを現金として保有し、償還請求に応じるか、それとも融資や株式、債券に投資するかを決定しなければならない。
小売業者(オンライン・実店舗を問わず)は、何十万もの製品の在庫を管理しなければならない。
自動車ディーラーは、顧客の需要を満たすために何台の車を保有すべきか決定しなければならない。
コンサルティング会社は、様々なコンサルティングプロジェクトの変動する需要を満たすために、何人の従業員を雇用しておくべきか決定しなければならない。

エネルギー貯蔵は、特に豊かな形態の在庫問題である。あり得るすべてのバリエーション(それは無数にある)を考察するわけではないが、本章の問題は以下の特徴を持つ。

送電網上の電力価格は非常に変動が大きい可能性がある。2000年代初頭には、典型的なエネルギー価格は1mwhあたり20~25ドル程度であったが、しばしば300ドルを超えるスパイクが生じ、極端な悪天候時には1,000ドルを超えることもあった。
風力エネルギーは予測可能であるが、あまり精度は高くない。これらの推定値を更新するためのローリング予測が利用可能である。
太陽光エネルギーは3種類の変動性を示す。すなわち、太陽が昇り沈むという日周サイクルの非常に予測可能な過程、非常に晴れた日や非常に曇った日の存在(これらは通常1日以上前から予測可能である)、そして1時間先でさえ予測が困難だが送電網に深刻な電力急増を引き起こしうる、局所的な雲の変動性である。
エネルギーの需要は変動するが、主に気温(そしてより小さい程度で湿度)に依存するため、比較的予測可能である。
エネルギーは現行の送電網価格で売買することができる。同様に、再生可能エネルギー源からのエネルギーは、(構成に応じて)現在の負荷(電力需要)を満たすために使用するか、貯蔵するか、あるいは送電網に売り戻すことができる。
送電網から到来する電力を交流(AC)から、電池に貯蔵するために必要な直流(DC)に変換する際には、概ね5~10パーセントの損失が生じる。

問題の枠組み化

3つの枠組み化に関する質問への答えは以下の通りである。

評価指標: 送電網から購入する電力に対して支払う期待価格を最小化したい。
決定: 各時間period(期間)ごとに送電網から購入する、あるいは送電網に売り戻すエネルギー量。
不確実性: 各時間期間における電力価格。

基本モデル

この問題については、様々なバリエーションを順に用いて、異なるモデリング上の課題を説明していく。まずは、電池を用いて価格の変動性を利用しながら送電網との間で売買を行う基本的なシステムから始める。この応用問題では、5分刻みの時間増分で時間を進めていく。これは、送電網上で価格が更新される頻度だからである(この時間刻み幅は送電網の運用者によって異なる)。

状態変数

本章の基本モデルでは、2つの変数のみを追跡すればよい。すなわち、時刻$t$において電池に貯蔵されているエネルギー量(メガワット時、すなわちmwh単位で計測)を表す$R_t$と、送電網上のエネルギー価格を表す$p_t$である。したがって、状態変数は次のようになる。

\[S_t = (R_t, p_t).\]

モデルに異なる要素を加えていくにつれ、状態変数はすぐにより複雑になっていく。例えば、価格を表す状態変数は、(以下で述べる遷移関数において)価格過程をどのようにモデル化するかに依存する。

決定変数

我々の唯一の決定は、送電網から購入するか、送電網に売却するかである。すなわち$x_t$は、送電網から購入する($x_t > 0$)、あるいは送電網に売却する($x_t < 0$)電力量である。

電池へエネルギーを移動させる、または電池からエネルギーを取り出す際には、その移動において得られるのは割合にして$\eta$のみであると仮定する。これは$1-\eta$の損失を意味する。簡単のため、この損失は電池を充電しているか放電しているかにかかわらず同じであると仮定する。

この決定は電池の容量によって制約されるため、次の制約条件を満たさなければならない。

\[x_t \leq \frac{1}{\eta} (R^{max} - R_t), \qquad x_t \geq -\eta R_t,\]

ここで、最初の制約は送電網から購入する場合($x_t > 0$)に適用され、2番目の制約は送電網に売却する場合($x_t < 0$)に適用される。

いつものように、決定は以下で定める方策$X^\pi(S_t)$によってなされると仮定する。

外生情報

本章の基本モデルにおいて、唯一の外生情報は価格の変化である。過去の価格に基づいて価格を予測するモデルなしに、各時間期間における価格が明らかになると仮定してよい。この場合、外生情報$W_t$は次のようになる。

\[W_{t+1} = p_{t+1}.\]

あるいは、価格の変化$\phat_t = p_t - p_{t-1}$を観測すると仮定することもでき、その場合は次のように書くことになる。

\[W_{t+1} = \phat_{t+1}.\]

遷移関数

状態変数の推移は次式で与えられる。

\[\begin{align} R_{t+1} &= \begin{cases} R_t + \eta x_t & x_t \geq 0, \\ R_t + \dfrac{x_t}{\eta} & x_t < 0. \end{cases} \label{eq:energytransition1}\\ p_{t+1} &= p_t + \phat_{t+1}. \label{eq:energytransition2} \end{align}\]

このように価格過程の変化を「観測する」というスタイルでモデル化することは、遷移関数を記述する際に有用である。実際には、通常は(変化ではなく)$p_{t+1}$を直接観測することになるため、その場合には明示的な遷移方程式は不要である。これら2つの式が、我々の遷移関数$S_{t+1} = S^M(S_t,x_t,W_{t+1})$を構成する。

後に、決定後状態$S^x_t$をモデル化することが有用であることが分かる。これは、決定$x_t$を下した直後だが、新しい情報が到着する前の状態である。決定後の資源状態は次のようになる。

\[R^x_t = \begin{cases} R_t + \eta x_t & x_t \geq 0, \\ R_t + \dfrac{x_t}{\eta} & x_t < 0. \end{cases}\]

貯蔵量を表す変数は決定論的に推移するため、次の決定前状態への遷移は単に次のようになる。

\[R_{t+1} = R^x_t.\]

一方、価格$p_t$は決定の影響を受けないため、決定後の価格は単に次のようになる。

\[p^x_t = p_t.\]

これにより、決定後状態は次のようになる。

\[S^x_t = (R^x_t, p_t).\]

目的関数

各期間において、得る(あるいは失う)金額は次式で与えられる。

\[C(S_t,x_t) = -p_t x_t.\]

したがって、目的関数は次式で与えられる標準的な問題となる。

\[\max_\pi \E \sum_{t=0}^T -p_t X^\pi(S_t),\]

ここで$S_{t+1} = S^M(S_t,X^\pi(S_t),W_{t+1})$は式$\eqref{eq:energytransition1}$および$\eqref{eq:energytransition2}$によって与えられる。さらに、初期状態$S_0$(すなわち$R_0$および$p_0$)を指定する必要があり、また$W_1, W_2, \ldots$を生成する方法が必要となる。これについて次に説明する。

不確実性のモデル化

第2章の資産売却問題では、価格を次式に従ってモデル化できると仮定した。

\[p_{t+1} = p_t + \varepsilon_{t+1},\]

ここで、$\varepsilon_{t+1}$は平均0、既知の分散を持つ正規分布に従うと仮定した。図8.3は、「地域限界価格」(電力業界の用語ではLMP)として知られる送電網の価格を1年間にわたって示したもので、送電網上の価格が示す途方もない変動性を示している。この変動性は、負荷の急増(あるいは電力の損失)によって短期的な不足が生じうることに起因する。需要は非弾力的である(送電網は負荷の100パーセントを満たすものとされている)ため、価格は短期間(価格は5分刻みで更新される)で20倍から50倍に跳ね上がることがある。

電力価格をモデル化する方法はいくつかある。以下では、この問題に用いられてきた4つの手法を説明する。

時系列モデル

時系列に関する文献は非常に豊富であるため、ここでは価格$p_{t+1}$を直近の価格履歴の関数として表す基本モデルのみを例示する。説明のため、直近3つの時間期間を用いることとし、モデルは次のように書ける。

\[\begin{align} p_{t+1} &= \thetabar_{t0} p_t + \thetabar_{t1} p_{t-1} + \thetabar_{t2} p_{t-2} + \varepsilon_{t+1}, \label{eq:energytimeseriesmodel}\\ &= \thetabar^T_t \phi_t + \varepsilon_{t+1}, \nonumber \end{align}\]

ここで

\[\phi_t = \begin{pmatrix} p_t \\ p_{t-1} \\ p_{t-2} \end{pmatrix}\]

は価格のベクトルである。ノイズ$\varepsilon \sim N(0,\sigma^2_\epsilon)$は、与えられた$\sigma^2_\epsilon$に対して次のように仮定する。

係数のベクトル$\thetabar_t = (\thetabar_{t0},\thetabar_{t1},\thetabar_{t2})^T$は再帰的に推定することができる。まず、係数ベクトルの初期推定値$\thetabar_0$から始めると仮定する。また、3行3列の行列$M_0$も必要となるが、今のところはスケーリングされた単位行列であると仮定してよい(より良いアイデアは以下で示す)。

$\thetabar_t$に対する基本的な更新式は次式で与えられる。

\[\thetabar_{t+1} = \thetabar_t - H_t\phi_t \varepsilon_{t+1},\]

誤差$\hat{\varepsilon}_t$は次式を用いて計算される。

\[\varepsilon_{t+1} = \thetabar^T_{t}\phi_t - p_{t+1}.\]

3行3列の行列$H_t$は次式を用いて計算される。

\[H_t=\frac{1}{\gamma_t}M_t,\]

行列 $M_t$ は次式を用いて再帰的に計算される。

\[M_t = M_{t-1} - \frac{1}{\gamma_t} (M_{t-1} \phi_t (\phi_t)^T M_{t-1}).\]

変数 $\gamma_t$ は次式を用いて計算されるスカラーである。

\[\gamma_t = 1 + (\phi_t)^TM_{t-1}\phi_t.\]

これらの式には $\thetabar_0$ と $M_0$ の初期推定値が必要である。これを行う一つの方法は、初期データを収集し、静的な推定問題を解くことである。$K$ 個の価格を観測すると仮定する。$Y_0$ を観測された価格 $p_3, p_4, \ldots, p_{K+3-1}$ からなる $K$ 要素の列ベクトルとする(モデルで直近3つの価格が必要となるため、3番目の価格から開始しなければならない)。

次に $X_0$ を $K$ 行からなる行列とし、各行が $p_k, p_{k-1}, p_{k-2}$ から構成されるものとする。$\thetabar$ の最良推定値は、正規方程式

\[\thetabar_0 = [(X_0)^T X_0]^{-1} (X_0)^T Y_0.\]

によって与えられる。最後に $M_0 = [(X_0)^T X_0]^{-1}$ とすると、行列 $M_t$ が時刻 $t$ における $[(X_t)^T X_t]^{-1}$ の推定値であることが分かる。

時間の経過に伴う変数間の関係を捉える時系列モデルの一群が存在する。もしこれらの手法を価格データにそのまま適用すると、結果はかなり悪いものとなるだろう。第一に、価格は正規分布に従わない。第二に、価格が負になることはあるものの、これはかなり稀である。しかし、このモデルを直接適用すると、分散 $\sigma^2_\epsilon$ がこの種のデータの高いノイズにキャリブレーションされていた場合、負の価格が生成される可能性が非常に高くなる。最後に、時間の経過に伴う価格の急上昇(ジャンプ)の挙動は、現実的なものとはならないだろう。

ジャンプ拡散

上記の線形モデルに対する主要な批判は、電力価格の研究でよく見られる大きなスパイクを捉えるのが不得手だという点である。この限界を克服するための単純なアイデアは、ジャンプ拡散モデルとして知られる手法を用いることである。これは、式 $\eqref{eq:energytimeseriesmodel}$ にもう一つのノイズ項を加えて次式を得るというものである。

\[p_{t+1} = \thetabar_{t0} p_t + \thetabar_{t1} p_{t-1} + \thetabar_{t2} p_{t-2} + \varepsilon_{t+1} + \mathbb{I}_t \varepsilon^J_{t+1}.\]

ここで、指示変数 $\mathbb{I}_t = 1$ はある確率 $p^{jump}$ で発生し、ノイズ $\varepsilon^J_{t+1}$ は平均 $\mu^{jump}$(典型的にはゼロよりかなり大きい)、分散 $(\sigma^{jump})^2$(かなり大きい)の正規分布に従う。

ジャンプ確率 $p^{jump}$、および平均と分散 $(\mu^{jump}, (\sigma^{jump})^2)$ を推定する必要がある。これは、$p^{jump} = 0$ となる基本モデルから始めることで行われる。この基本モデルを用いて $\sigma^2_\epsilon$ を推定する。次に、標準偏差3つ分といった許容範囲(すなわち $3 \sigma_\epsilon$)を選び、この範囲外にある観測値は別のノイズ源によるものとみなす。$p^{jump}$ を、このような観測値が発生する時間区間の割合とする。そして、これらの観測値の平均と標準偏差を計算して $(\mu^{jump}, (\sigma^{jump})^2)$ を得る。

ここで終わってはいけない。データからこれらの極端な変動を取り除いた後、これらの観測値を除いて線形モデルを再び当てはめるべきである。標準的な手法は、これらの推定値が変化しなくなるまでこの処理を数回繰り返すことである。

ジャンプ拡散モデルは裾野をより良く再現するが、それでも正規分布の裾野の挙動に依存している。より良い当てはめは、ノイズの分散が気温、特に極端な気温に依存していることを認識することで得られる。気温を、氷点下、華氏90度以上、その中間の3つの範囲にグループ化することが考えられる。気温依存性を導入することは、状態変数の集合に別の変数を追加することになり、追加の複雑さをもたらす(この影響は方策のクラスに依存する)。

分位分布

他のパラメトリック分布を当てはめることも可能かもしれないが、強力な戦略は、データから累積分布を数値的に計算し、いわゆる分位分布を作成することである。これを計算するには、単に価格を小さい順から大きい順に並べ替える。この順序付けられた系列を $\ptilde_t$ で表す。ただし $\ptilde_{t-1} \leq \ptilde_t$ とする。$T = 105,210$ を年間の5分間隔の時間区間数とする。価格が $\ptilde_t$ より小さい時間区間の割合は $t/T$ となる。次式を用いて累積分布を作成できる。

\[F_P(\ptilde_t) = \frac{t}{T},\]

これは図 8.4 に示されている。

任意の $p$ に対して連続分布 $F_P(p)$ を作成するには、最大の $\ptilde_t < p$ を見つけ、$F_P(p)$ をこの値に等しく設定することで階段関数を作ることができる。関数 $F_P(p)$ は、既知のパラメトリックな形式に分布を当てはめているわけではないため、ノンパラメトリック分布の一形態である。良い点は、これがデータに完全に一致するため、電力価格に生じる極端な裾野を正確に表現できることである。欠点は、これらの分布を作成するには良質なデータセットが必要であり、分布をパラメトリックモデルで当てはめる場合のように少数のパラメータだけを保存するのではなく、データセット自体を保持して分布を計算しなければならないことである。

この分布からサンプリングするには、0から1の間で一様分布する乱数変数 $U$ を生成すればよい。例えば $U= 0.70$ を生成したとしよう。すると、図 8.4 に示されているように、$F_P(p^{.70}) = 0.70$ に対応する価格 $p^{.70}$ を見つけたい。これを数学的に表現するには、$F_P(p) = u$ を生成する価格 $p$ を返す逆関数 $F^{-1}_P(u)$ を定義する。一様乱数変数 $U$ を単にサンプリングし、価格 $p=F^{-1}_P(U)$ を観測することで、価格分布から繰り返しサンプリングできる。

変換データを用いたハイブリッド時系列

強力な戦略は、経験分布の使用を古典的な時系列手法と組み合わせることである。まず、価格データに経験分布を当てはめ、累積分布 $F_P(p)$ を得る。次に、$p_t$ を価格とし、$u_t = F_P(p_t)$ を計算する。ただし $0\leq u_t \leq 1$ である。次に、$\Phi(z)$ を標準正規確率変数 $Z \sim N(0,1)$ の累積分布とし、$\Phi^{-1}(u)$ をその逆関数とする。次に $z_t = \Phi^{-1}(u_t)$ とする。$p_t \rightarrow u_t \rightarrow z_t$ の写像プロセスは図 8.5 に示されている。

この方法を用いて、非正規性の強い価格 $p_t$ を、平均0、分散1の正規分布に従う値の系列 $z_t$ に変換することができ、相関も捉えることができる。そして、この系列 $z_t$ に対して任意の時系列モデリングを行うことができる。その後、この正規化されたモデルから得られる推定値は、図 8.5 の経路を逆方向にたどることで価格に戻すことができる。すなわち $u_t = \Phi(z_t)$、そして $p_t = F^{-1}_P(u_t)$ である。

この戦略は、正規分布に従わないデータを扱う際に非常に有効であり、金融分野で人気のあるジャンプ拡散モデルよりもはるかに優れた性能を発揮する。

方策の設計

この問題を解くために、2種類の方策のクラスとハイブリッドを一つ用いて説明する。

方策探索 – バイローセルハイ(安く買って高く売る)というパラメータ化された単純な方策を使用する。これはPFA(policy function approximations、方策関数近似)のクラスに属する方策である。
先読み方策 – ベルマン方程式を用いて、決定を今下すことが将来に与える影響を近似する価値関数近似を生成する。これはVFA(価値関数近似に基づく方策)のクラスに属する。
ハイブリッド方策 – 最後に、価値関数に基づく方策から始まり、方策探索へと切り替えてさらに方策を調整する、調整可能な方策のクラスを導入する。これは最初はVFA方策として始まるが、方策探索を用いて価値関数近似として始まったもののパラメータを調整する際にパラメトリックなコスト関数近似(CFA)へと移行する。

ベルマン方程式に基づく方策では、状態 $S_t$ にあり、決定 $x_t$ を下し、その後ランダムな外生情報 $W_{t+1}$ を観測することから生じる値 $V_{t+1}(S_{t+1})$ を計算(または近似)する必要がある。これらの手法は最短経路問題の文脈で最初に見た。今回の大きな違いの一つは、状態 $S_{t+1}$ が $S_t$ と $x_t$ を条件としてランダムであることである(最短経路問題では、コスト $\chat_t$ のみがランダムであった)。また、今回の状態変数は、単なる離散的なノードではなく、2つの連続次元を持っている。

まずバイローセルハイ方策を説明し、その後ベルマン方程式の近似に基づく3つの手法を紹介する。

最適な方策を生成する古典的な後向き動的計画法。
後向き近似動的計画法。
前向き近似動的計画法。

最後に、ベルマン方程式からの価値関数近似と方策探索の一形態を組み合わせたハイブリッド方策の説明で締めくくる。

バイローセルハイ

バイローセルハイ方策は、価格が下限を下回ったときに電池を充電し、上限を上回ったときに売却するという単純な原則で動作する。この方策は次のように書ける。

\[X^{low-high}(S_t\vert \theta) = \begin{cases} -1 & \text{if } p_t \leq \theta^{buy}, \\ 0 & \text{if } \theta^{buy} < p_t < \theta^{sell}, \\ +1 & \text{if } p_t \geq \theta^{sell}. \end{cases}\]

次に $\theta = (\theta^{buy}, \theta^{sell})$ を調整する必要がある。価格 $p_t(\omega)$ のサンプルパスに従って方策を評価する(あるいは価格変化 $\phat(\omega)$ を観測している場合もある)。数学モデルからサンプルパスを生成すると仮定すると、サンプルパス $\omega^1, \ldots, \omega^N$ を生成できる。そして、各サンプルパスにわたって方策の性能をシミュレートし、次式を用いて平均を取ることができる。

\[\Fbar^{low-high} = \frac{1}{N} \sum_{n=1}^N C\big(S_t(\omega^n),X^{low-high}(S_t(\omega^n)\vert \theta)\big).\]

$\theta$ を調整するには、次の問題を解く必要がある。

\[\begin{align} \max_\theta \Fbar^{low-high}(\theta\vert S_0). \label{eq:buylowpolicysearch} \end{align}\]

$\theta$ は2次元しか持たないため、一つの戦略は各次元を離散化し、2つの次元のすべての可能な値を探索する完全グリッド探索を行うことである。一般的な離散化は、範囲を5パーセント刻みに分割することである。境界を含めると、各パラメータの21個の値を表現する必要があり、これは大きさ441点のグリッドを作成することになるが、多くの問題では(自明ではないものの)管理可能である。

ブルートフォースのグリッド探索は、推定値 $\Fbar^\pi(\theta)$ の分散が比較的小さくなるように十分な回数のシミュレーション $N$ を実行した場合にのみ機能することに注意する。しかし、第7章で紹介したように、方策の性能に関するノイズの多い推定値であっても $\theta$ の探索を行う手法が存在する。

最も単純な方策には常に調整可能なパラメータが特徴として存在するため、$\eqref{eq:buylowpolicysearch}$ で与えられる最適化問題は繰り返し現れることになる。$\theta$ に関する $\Fbar^{low-high}(\theta\vert S_0)$ の導関数を計算(または近似)できる場合、問題 $\eqref{eq:buylowpolicysearch}$ は導関数に基づく手法を用いて解くことができる。これが不可能な場合は、導関数を用いない手法を使わなければならないが、これはまさに第2章で直面した問題である。

後向き動的計画法

後向き動的計画法は、ベルマン方程式

\[\begin{align} V_t(s_t) = \max_{x_t} \left(C_t(s,x_t)+ \E\{V_{t+1}(S_{t+1})\vert S_t,x_t\} \right), \label{eq:energystoragebellman} \end{align}\]

を直接解くことを伴う。ここで $S_{t+1} = S^M(s_t,x_t,W_{t+1})$ であり、期待値は確率変数 $W_{t+1}$ に関して取られる。$W_{t+1}$ が離散的であり、$\Wcal = \lbrace w_1, w_2, \ldots, W_M\rbrace $ の中の値を取ると仮定し、確率分布を次式を用いて表現する。

\[f^W(w\vert s_t,x_t) = Prob[W_{t+1} = w\vert s_t,x_t].\]

この分布を状態 $s_t$ と決定 $x_t$ に依存するものとして記述しているが、これは問題によって異なる。例えば、価格の変化 $p_{t+1}-p_t$ が現在の価格 $p_t$ に依存する(価格が非常に高ければ下落する可能性が高い)と合理的に仮定できるかもしれず、これは $s_t$ に条件付けする理由となる。電力網から大量の電力を購入することが価格を押し上げる場合、$x_t$ への依存性さえ必要になるかもしれない。

これにより、式 $\eqref{eq:energystoragebellman}$ を次のように書き直すことができる。

\[V_t(s_t) = \max_{x_t} \left(C_t(s_t,x_t)+ \sum_{w\in\Wcal} f^W(w\vert s_t,x_t)V_{t+1}(S^M(s_t,x_t,w)) \right).\]

後向き動的計画法の素朴な実装には、4つのループが現れる。

時刻 $T$ から時刻 $0$ まで時間を後向きにステップするループ。
可能なすべての状態 $s_t\in\Scal$ にわたるループ(より正確には、これは時刻 $t$ における状態 $S_t$ の可能な値の集合である)。
最大化問題を解くために可能なすべての決定 $x_t$ を探索する際に必要となるループ。
$V_t(s)$ を計算するために必要な総和に含まれる、確率変数 $W$ の可能なすべての値にわたるループ。

後向き動的計画法

ステップ0. 初期化： すべての状態$S_{t+1}$について、終端貢献$V_{T+1}(S_{T+1})=0$を初期化する。

ステップ1. $t=T, T-1, \ldots, 1, 0$について実行：

ステップ2. すべての$s\in\Scal$について、以下を計算：

$$V_t(s_t) = \max_{x_t} \left(C_t(s_t,x_t)+ \sum_{w\in\Wcal} f^W(w\vert s_t,x_t)V_{t+1}(S^M(s_t,x_t,w)) \right)$$

各ループが取り得る値の範囲を考えてみることは有用である。エネルギー問題では、1日を1時間刻みで蓄電デバイスを最適化するとすれば、24の時間ステップが得られる。5分刻みの時間ステップを用いる場合（一部の系統運用者は5分ごとに価格を更新する）、24時間の計画期間は288の時間期間を意味する（1週間にわたって計画する場合はこれを7倍する）。周波数調整を行う場合は、2秒ごとに決定を下す必要があり、これは1日あたり43,200の時間期間に相当する。

我々の状態変数は$S_t = (R_t,p_t)$から構成されるため、すべての状態にわたるループを、$R_t$のすべての値、そして$p_t$のすべての値にわたる入れ子のループに置き換える必要がある。両者は連続であるため、それぞれ離散化しなければならない。資源変数$R_t$は、1つの時間増分でどれだけ充電または放電し得るかに基づいて増分に分割する必要がある。次に系統価格$p_t$を離散化しなければならない。系統価格は（極端な場合）低くはマイナス100ドル、高くは10,000ドルにまでなり得る。妥当な戦略としては、経験分布を構築し、累積分布の例えば2パーセントごとの増分に対応する価格を表現することが考えられ、これにより50通りの価格が得られる。

充電・放電決定の数は3つ（充電、放電、何もしない）と少ない場合もあれば、異なる速度で充電・放電できる場合ははるかに多くなり得る。

最後に、確率分布$f^W(w)$は価格のランダムな変化$\phat_{t+1}$の分布である。ここでも、$\phat_{t+1}$の変化の経験分布を構築し、累積分布を例えば2パーセントの増分に離散化することを推奨する。

我々の基本モデルのように2次元の状態変数を持つ場合、既に5つのループがある（時間、2つの状態変数、$x$に対するmax演算子、そして$W$の結果にわたる総和）。これは計算コストが高くなり、まだ始まったばかりである。ここで、式$\eqref{eq:energytimeseriesmodel}$の時系列モデルを使用しており、価格$(p_t, p_{t-1}, p_{t-2})$を追跡する必要がある場合を想像してほしい。この場合、我々の状態変数は

\[S_t = (R_t, p_t, p_{t-1}, p_{t-2}).\]

となる。この場合、7つの入れ子ループを持つことになる。複雑さは連続変数の離散化に依存するが、この後向き動的計画法アルゴリズムを実行するには1年（あるいはそれ以上）を要する可能性がある。

この比較的単純な問題に対してさえベルマン方程式を用いることがいかに困難であるかを考えると、この特定のアプローチが今なお授業で教えられていることは驚くべきことである。この複雑さを踏まえ、ベルマン方程式を近似する手法について広範な研究が行われてきており、これらは近似動的計画法や強化学習といった名称の下にまとめられている。ここでは、後向きADPと前向きADPと呼ばれる、ベルマン方程式を近似する2つの戦略を説明する。

後向き近似動的計画法

強力なアルゴリズム戦略として「後向き近似動的計画法」として知られるものがある。このアプローチは、1つの違いを除いて、上記でちょうど行ったのと同じように進行する。すべての状態にわたってループする代わりに、ランダムなサンプル$\Shat$を選択する。次に、当初行ったのとまったく同様に状態$s\in\Shat$にある価値を計算し、対応する値$\vhat$を計算する。これを$N$回繰り返し、データセット$(\shat^n, \vhat^n), n=1, \ldots, N$を取得すると仮定する。次に、これを用いて、例えば次式で与えられる線形モデルのような統計モデルを当てはめる：

\[\begin{align} \Vbar(s) = \theta_0 + \theta_1 \phi_1(s) + \theta_2 \phi_2(s) + \ldots + \theta_F \phi_F(s), \label{eq:energylinearvfa} \end{align}\]

ここで、$\phi_f(s), f=1, \ldots, F$は適切に選択された特徴量の集合である。特徴量の例としては

\[\begin{align*} \phi_1(s) &= R_t, \\ \phi_2(s) &= R^2_t, \\ \phi_3(s) &= p_t, \\ \phi_4(s) &= p^2_t, \\ \phi_5(s) &= p_{t-1}, \\ \phi_6(s) &= p_{t-2}, \\ \phi_7(s) &= R_t p_t. \end{align*}\]

が挙げられる。定数項$\theta_0$を含めても、このモデルは推定すべき係数がわずか8個である。数百の状態をサンプリングすれば、良好な統計的近似を得るには十分すぎるほどである。この方法論は状態変数の数に対して比較的鈍感であり、もちろん変数のいずれかが連続であっても問題はない。

上記の線形モデルのようなパラメトリックモデルを使用する際に常に課題となるのは、特徴量$\phi_f(S_t)$を指定しなければならないことである。ニューラルネットワークが普及するにつれ、研究者たちはこのアプローチを、数百万のパラメータを推定する必要があるかもしれない深層ニューラルネットワークを含めて用いるようになった。このアプローチの利点は、モデルの構造を指定する必要がなくなることであるが、その代償としてはるかに多くの観測値が必要になる。深層ニューラルネットワークは任意の関数を近似できるという魅力的な性質を持つが、これはノイズをモデル化してしまう可能性があることも意味する。ニューラルネットワークはまた、単調性（在庫が大きいほど価値が大きい）や凸性といった既知の問題構造を再現することにも苦労する。

我々は、後向きADPが少数の問題群において非常によく機能することを見出してきた（後向きADPとベンチマークとの比較の要約については、Reinforcement Learning and Stochastic Optimizationの15.4節を参照）が、保証は一切なく、その性能は明らかに効果的な特徴量の集合を選択できるかどうかに依存する。あるアプリケーションでは、標準的な後向きMDPアルゴリズムの実行時間30日を20分にまで短縮し、（1か月の実行時間で生成された）最適解の5パーセント以内の解を得た。しかし、繰り返すが、この性能が保証されるわけではない。

前向き近似動的計画法

前向き近似動的計画法は直観的な方法で機能する。アルゴリズムの最初の$n-1$回の反復から計算した決定後状態$S^x_t$周りの価値関数近似$\Vbar^{x,n-1}_t(S^x_t)$から始めると想像してほしい。決定後状態という考え方は第1章で最初に紹介したが、これは決定を下した直後、しかし新しい情報が到着する前の状態である。

さて、アルゴリズムの$n$回目の反復において、時間を前進させる際のサンプリングを導くサンプルパス$\omega^n$に従い、特定の状態$S^n_t$にいると想像してほしい。決定後状態へと導く関数$S^{x,n}_t = S^{M,x}(S^n_t,x)$があると仮定する。$S^n_t = (R^n_t,p^n_t)$である我々のエネルギー問題については、決定後状態は

\[S^{x,n} = (R^n_t+x^n_t, p^n_t).\]

となる。次に、

\[x^n_t = \argmax_x \big(C(S^n_t,x) + \Vbar^{x,n-1}_t(S^{x,n}) \big).\]

を用いて決定を下す。$S^n_t$と我々の決定$x^n_t$が与えられれば、$W_{t+1}(\omega^n)$をサンプリングし、これは価格の変化$\phat^n_{t+1}$に変換される。次に、次の状態へとシミュレートしていく

\[S^n_{t+1} = (R^n_t+x^n_t, p^n_t + \phat^n_{t+1}(\omega)).\]

このように、我々は時間を前方にシミュレートしているだけなので、状態変数がどれほど複雑であるかは気にする必要がない。価値関数近似$\Vbar^{n-1}_t$を更新するにはさまざまな戦略がある：

前向き近似動的計画法

ステップ0. 初期化： $V^{\pi,0}_t,~t\in\Tcal$を初期化する。$n = 1$を設定する。$S^1_0$を初期化する。

ステップ1. $n = 1, 2, \ldots, N$について実行：

ステップ2. $m = 1, 2, \ldots, M$について実行：

ステップ3. サンプルパス$\omega^m$を選択する。

ステップ4. $\vhat^m = 0$を初期化する。

ステップ5. $t = 0, 1, \ldots, T$について実行：

ステップ5a. 以下を解く：

$$x^{n,m}_t = \argmax_{x_t\in\Xcal^{n,m}_t} \big(C_t(S^{n,m}_t,x_t) + V^{\pi,n-1}_t(S^{M,x}(S^{n,m}_t,x_t))\big)$$

ステップ5b. 以下を計算：

$$S^{x,n,m}_t = S^{M,x}(S^{n,m}_t,x^{n,m}_t), \qquad S^{n,m}_{t+1} = S^M(S^{x,n,m}_t,x^{n,m},W_{t+1}(\omega^m)).$$

ステップ6. $t = T-1,\ldots, 0$について実行：

ステップ6a. 経路コストを（$\vhat^m_{T} = 0$を用いて）累積する：

$$\vhat^m_t = C_t(S^{n,m}_t,x^m_t) + \vhat^m_{t+1}$$

ステップ6b. 時間$t$から始まる方策の近似価値を更新する：

$$\Vbar^{n,m}_{t-1} \leftarrow U^V(\Vbar^{n,m-1}_{t-1}, S^{x,n,m}_{t-1}, \vhat^m_t)$$

ここで、通常は$\step_{m-1} = 1/m$を用いる。

ステップ7. すべての$t = 0, 1, \ldots, T$について、方策の価値関数$V^{\pi,n}_t(S^x_t) = \Vbar^{n,M}_t(S^x_t)$を更新する。

ステップ8. 価値関数$(V^{\pi,N}_t)_{t=1}^T$を返す。

これは、価値関数をどのように近似するかに依存するため、実際の更新を更新関数$U^V(\cdot)$に委ねている。

前向き近似動的計画法は、高次元の問題にスケールするという点で魅力的である。すべての可能な状態や結果にわたってループすることは一度もない。実際、価値関数を適切に近似して強力なアルゴリズムを活用できるようにすれば、高次元の決定$x$すら扱うことができる。しかし、前向きADPは（後向きADPと同様に）性能保証をほとんど持たない。

方策探索とVFA方策のハイブリッド

価値関数をどのように近似することを選択しても、我々の方策は次式で与えられる

\[X^{VFA}(S_t) = \argmax_x \big(C(S_t,x) + \Vbar^x_t(S^x_t)\big).\]

方策を前向きにシミュレートする際には、$\omega^n$に外生情報（すなわち、価格の変化の集合）のサンプルパスを表させる。しばしば、我々は履歴データに対して方策をテストすることになり、その場合サンプルパスは1つだけである。しかし、不確実な価格に関する数理モデルを開発している場合には、方策の価値を近似するために用いるサンプルパス$\omega$を生成することができる（複数のサンプルパスを生成して平均を取ることも可能である）：

\[\Fbar^{VFA}(\omega\vert S_0) = \sum_{t=0}^T C\big(S_t(\omega),X^{VFA}(S_t(\omega))\big).\]

これは、この方法がロジスティクスにおいて生じ得るような高次元の問題にすら適合していることを意味する。

価値関数近似を構築している場合（これは先読み方策のクラスに属する）、通常、方策を調整するステップはもはや存在しない。しかし、これは試みることができないという意味ではない。我々の価値関数が式$\eqref{eq:energylinearvfa}$の線形モデルで与えられると仮定する。このとき、我々の方策は次式で書くことができる

\[X^{VFA}(S_t\vert \theta) = \argmax_x \left(C(S_t,x) + \sum_{f=1}^F \theta_f \phi_f(S_t)\right).\]

我々の後向きまたは前向きADPアルゴリズムの1つを用いて$\theta$の初期推定値を得るのは理にかなっているが、上述したように、結果として得られる方策が高品質であるという保証はない。しかし、これを出発点として用いることで、常により良いものにすることができる、

\[\Fhat^{VFA}(\theta,\omega\vert S_0) = \sum_{t=0}^T C\big(S_t(\omega),X^{VFA}(S_t(\omega)\vert \theta)\big).\]

さて、次のように定式化できる方策探索問題を解くだけでよい

\[\begin{align} \max_\theta \Fbar^{VFA}(\theta,\omega\vert S_0). \label{eq:optthetavfa} \end{align}\]

方策探索の初期の例では、$\theta$はスカラーであったため、この問題は比較的容易であった。ここでは、$\theta$は数十次元を持ち得るベクトルである。この問題については後で改めて取り上げる。

ADPに関するいくつかの注意点

我々はこの問題設定を用いて、状態にある価値を近似するという考え方に基づく手法について比較的深く掘り下げたツアーを提供してきた。これは「近似動的計画法」あるいは「強化学習」といった用語のもとで研究されてきた。これらの手法は学術研究コミュニティから多大な注目を集めてきたが、実務においては、これらの手法は容易ではない。逐次決定問題はいたるところに存在するが、実務における成功事例は比較的稀である。

各離散（または離散化された）状態の価値を推定するルックアップテーブル近似は、状態変数が3次元を超えると拡張性がない。我々の線形近似のような近似戦略を用いることは、通常うまく機能しない。なぜなら、これらの近似は、任意の状態を訪れる可能性があるため、大域的に正確でなければならないからである。同時に、局所近似（これはノンパラメトリックモデルの一形態である）は、局所近似の柔軟性が不安定性をもたらすため、苦戦することがある。

このプロセスをさらに複雑にしているのは、決定を下すために近似価値関数に依存しているという点であり、これが悪循環を生み出す。我々の初期の近似はあまり良くなく、その結果、悪い決定につながる。これらの悪い決定が今度は価値関数近似の更新に用いられ、そこから下降スパイラルが見えてくる。

価値関数近似をチューニングするというアイデアは、式$\eqref{eq:optthetavfa}$で行ったように、方策の性能を直接最適化するため有望である。奇妙なことに、このアイデアは広く使われていない。ここでは、これらのシミュレーションを用いて$\theta$をチューニングすることは容易ではないということだけを指摘しておく。したがって、これらのアプローチを試そうとする読者には注意を促しておきたい。

何を学んだか

第1章からの単純な在庫問題を、物理的変数と情報的変数の両方を伴うエネルギー貯蔵の文脈で再訪した。
不確実性モデリングの豊かさを示すために、電力価格に関するさまざまな確率モデルを説明した。
PFA(安値で買い、高値で売る)、VFA方策(後ろ向き近似動的計画法に基づく)、および前向き近似動的計画法を含む、さまざまな方策を説明した。
最後に、近似動的計画法の手法を用いて価値関数近似の線形モデルを推定し、その線形モデルの係数に対して直接方策探索を行うというアイデアを導入した。したがって、これは当初はVFAベースの方策であるが、パラメータ化された目的関数を持つCFA方策の一形態へと転じることになる。

演習問題

復習問題

電力価格が裾の重い(ヘビーテイルな)性質を持つことを踏まえると、時系列モデルを用いることの何が問題なのか。
通常の変動範囲内(標準偏差3つ分)に収まる価格と、より極端な観測値とをどのように分離するか、簡潔に概説せよ。
過去のデータを用いて経験分布をフィッティングすれば、履歴に一致する確率分布が得られるはずである。それでもなお、価格の確率モデルにはどのような誤差が残り得るか。
ハイブリッド時系列の節におけるデータ変換が何を達成しているのか、言葉で説明せよ。
古典的な後ろ向き動的計画法は、次元の呪いによってすぐに破綻する。後ろ向き近似動的計画法がどのように次元の呪いを克服するのかを言葉で説明せよ。例えば、状態変数の次元数を2倍にした場合、これが後ろ向きADPをどのように複雑にするかを説明せよ。

問題解決演習

本文で示されている通りに、エネルギー貯蔵問題の基本モデルの5つの要素を書き出せ。方策が安値で買い、高値で売る方策であると仮定して、目的関数を書き出せ。 $$ X^{low-high}(S_t\vert \theta) = \begin{cases} -1 & \text{if } p_t < \theta^{buy}, \\ 0 & \text{if } \theta^{buy} \leq p_t \leq \theta^{sell}, \\ +1 & \text{if } p_t > \theta^{sell}. \end{cases} $$ 方策のパラメータについて探索を行う形で目的関数を書き出せ。また、目的関数中の期待値を、各確率変数を反映した入れ子形式を用いて書き出せ。
不確実性のモデリングに関する節では、以下で与えられる価格の時系列モデルを導入した。 $$ p_{t+1} = \thetabar_{t0} p_t + \thetabar_{t1} p_{t-1} + \thetabar_{t2} p_{t-2} + \varepsilon_{t+1}. $$ 本書では係数ベクトル$\thetabar_t = (\thetabar_{t0},\thetabar_{t1}, \thetabar_{t2})$の更新式を説明している。状態$S_t$が時刻$t$以降のシステムをモデル化するために必要な*すべて*の情報から成ることを踏まえ、この価格プロセスを扱うための更新された状態変数と遷移関数を示せ。

プログラミング演習

これらの演習では、tinyurl.com/sdagithubにあるPythonモジュールEnergyStorage_Iを使用する。

このPythonモジュールを用いて、パラメータベクトル$\theta = (\theta^{buy}, \theta^{sell})$に対するグリッドサーチを実行せよ。$\theta^{sell}$を価格に対して1.0ドルから100.0ドルまで＄1刻みで変化させ、$\theta^{buy}$を1.0ドルから$\theta^{sell}$まで同様に＄1刻みで変化させよ。価格は8日間の実際の過去の時間毎の価格を用いる。
上述の後ろ向き動的計画法の戦略(このアルゴリズムはすでにPythonモジュールに実装されている)を用いて、最適方策を求めよ。価格プロセスは以下に従って推移すると仮定する。 $$ p_{t+1} = p_t + \varepsilon_{t+1}, $$ ここで$\varepsilon_{t+1}$は、実際の過去の価格の差分に基づく経験分布に従う。
1. 価格を＄1刻み、次に＄0.50刻み、最後に＄0.25刻みで離散化してアルゴリズムを実行せよ。3つの離散化水準それぞれについて状態空間のサイズを計算し、状態空間のサイズに対して実行時間をプロットせよ。
2. ＄1の離散化における最適価値関数を用いて、(a)で見つけた最良の買い売り方策との性能を比較せよ。
[tinyurl.com/sdamodelingsupplements](https://tinyurl.com/sdamodelingsupplements/)からスプレッドシート「Chapter8_electricity_prices」をダウンロードせよ。以下の問いには「electricity prices」タブのデータを用いること。
1. データセット中の1週間にわたる特定の時間帯についてランダムに選ばれた価格を$P$として、経験累積分布$F_P(p) = Prob[P \leq p]$を構築せよ。
2. $F^{-1}_P(u)$を逆累積分布とし、$u$は0から1の間の値とする。$u = 0, 0.1, 0.2, \ldots, 0.9, 1.0$に対応する価格$p(u) = F^{-1}_P(u)$を求めよ。これらの各価格に1/11の確率を与えて累積分布を求め、(a)で作成した累積分布と比較せよ。両者は一致するように見えるか。
スプレッドシート「Chapter8_electricity_prices」を用いて、以下の形式の平均回帰モデルをフィッティングせよ。 $$ \begin{align} p_{t+1} = p_t + \beta (\mubar_t - p_t) + \varepsilon_{t+1} \label{eq:priceexercise} \end{align} $$ ここで $$ \mubar_t = (1-\alpha)\mubar_{t-1} + \alpha p_t. $$ $\alpha = 0.15$と仮定する。以下を最小化する$\beta$を求めよ。 $$ G(\beta) = \sum_{t=0}^T \big(p_{t+1} - (p_t + \beta (\mubar^t-p_t))\big)^2. $$
1. 単純な一次元探索(例えば0から1まで0.1刻みで値を試す)によって平均回帰モデルをフィッティングせよ。
2. 標本から$\varepsilon$の標準偏差$\sigma$を計算せよ(平均は0であると仮定する)。ここでは、単一の一定した標準偏差を仮定する一方で、平均$\mubar_t$は時間とともに変化することを許容している点に注意せよ。
3. (a)で見つけた$\beta$の値を用い、平均0・標準偏差$\sigma$の正規分布から$\varepsilon_{t+1}$をサンプリングすることで、式$\eqref{eq:priceexercise}$を用いてサンプルパスを10本生成せよ。サンプルパスをグラフにプロットし、そのふるまいを過去の価格と比較せよ。似ているように見えるか。
次に、以下で与えられるジャンプ拡散モデルをフィッティングする。 $$ p_{t+1} = p_t + \beta(\mubar_t - p_t) + \varepsilon_{t+1} + J_{t+1} \varepsilon_{t+1}, $$ ここで$J_{t+1} = 1$は(以下で計算する)あるジャンプ確率のもとで発生し、それ以外は0であり、$\varepsilon^J_{t+1}$はジャンプが発生した際のジャンプ幅を表すランダムな量である。ジャンプ拡散モデルをフィッティングし、その結果を履歴と比較するために、以下の手順に従うこと。
1. 演習問題11で得た$\beta$の値を用い、データを走査して範囲$[\mubar_t \pm 3 \sigma]$の外側に外れるすべてのデータ点を特定せよ。
2. 演習問題11で見つけたのと同じ$\beta$の値を用いて、データに対して再度平均回帰を実行せよ。ただし今回は、(a)で除外されなかったデータ点のみを含めること。除外されなかったデータについて、$\sigma$の新たな平均と標準偏差を求めよ。
3. 保持されたデータ点に対して、(b)をもう一度繰り返し、再び$\pm 3 \sigma$の範囲外にあるデータ点を除外せよ。
4. (c)を終えた時点で除外された点の割合として、ジャンプの確率を計算せよ(この時点でデータ点を除外する処理を2回実行したことになる)。また、除外された点の平均と標準偏差も計算せよ。
5. 保持された点と除外された点それぞれについて求めた最終的な平均と分散を用い、かつジャンプ拡散確率を用いてジャンプをサンプリングしながら、ジャンプ拡散モデルのシミュレーションを10回実行せよ。これらのシミュレーションを履歴と比較し、得られた価格パスが演習問題11で得たものよりも現実的であるかどうかを論じ、価格パスを実際の履歴と比較せよ。
変換された価格を用いて演習問題11および12の一部を繰り返すことで、価格の良好なフィッティングを再度試みる。
1. 演習問題10の累積分布を用い、恒等式$U_t = F_P(p_t)$を用いて各価格を一様分布に従う確率変数に変換せよ。
2. 次に、一様分布に従う確率変数$U_t$を、$Z_t = \Phi^{-1}(U_t)$を用いて平均0・分散1の正規分布に従う確率変数に変換せよ。ここで$\Phi(z)$は標準正規分布の累積分布であり、$\Phi^{-1}(U_t)$はその逆関数である。これはExcelのnorm.s.inv(p)関数によって捉えられ、確率$p$(変数$U_t$で与えられる)に対応する$Z$値を返す。
3. ここで、演習問題11の平均回帰方程式における$\beta$を再びフィッティングする必要がある。ただし今回は、価格$p_t$の代わりに正規化された量$Z_t$を用いる。それ以外はすべて同じである(したがって、生の価格に対して$\beta$をフィッティングした際の手順にそのまま従えばよい)。
4. (新たな$\beta$の値を持つ)(c)のモデルを用いて、$Z_t$値のサンプルパスを作成せよ。次に、$U_t = \text{norm.s.dist}(Z_t,1)$を用いて$Z_t \leq z$(0から1の間で一様分布する)となる確率を求めよ。最後に、演習問題11で求めた累積分布を用いて、$U_t$値を価格に逆写像せよ。1本のサンプルパスをプロットせよ(Excelに習熟していればそれほど難しくないが、そうでなければこの最後のステップが厄介な部分である)。
5. 得られたサンプルパスのふるまいを、過去の分布と比較せよ。価格の分布は完璧に一致するはずであるが、価格の系列自体は依然として良好なフィッティングには見えないかもしれない点に注意せよ。それでもなお、どのような誤りを犯している可能性があるか。