ファインマン・カッツの公式による確率過程のシミュレーション

この例では、資産の予想最終価格を表す偏微分方程式を求めます。この価格は、確率微分方程式によって与えられる確率過程です。

以下の手順が含まれます。

確率微分方程式を使用したモデルのパラメーターの定義
伊藤ルールの適用
ファインマン・カッツ方程式の求解
期待される資産売却の時間の計算

1.確率微分方程式を使用したモデルのパラメーターの定義

時間間隔[0,T]で定義される資産価格のモデルX(t)は、 $d X = μ (t, X) d t + σ (t, X) d B (t)$ という形式の確率微分方程式によって定義される確率過程です。ここで、B(t)は単位分散パラメーターを持つウィーナー過程です。

シンボリック変数Tと次のシンボリック関数を作成します。

r(t)はスポット金利を表す連続関数です。このレートによって、時間Tにおける最終的な支払いの割引係数が決定します。
u(t,x)は、割引された先物価格の期待値です。これは、X(t) = xの条件下で $X (T) \exp (- \int_{t}^{T} r (t) d t)$ として計算されます。
$μ (t, X)$ と $σ (t, X)$ は、確率過程X(t)のドリフトと拡散です。

symsmu(t, X)sigma(t, X)r(t, X)u(t, X)T

ファインマン・カッツの定理によると,uは偏微分方程式 $\frac{\partial u}{\partial t} + \frac{σ^{2}}{2} \frac{\partial^{2} u}{\partial x^{2}} + μ \frac{\partial u}{\partial x} - u r = 0$ を時間Tにおける最終条件で満たします。

eq = diff(u, t) + sigma^2*diff(u, X, X)/2 + mu*diff(u, X) - u*r;

数値ソルバーpdepeは初期条件を処理します。最終条件を初期条件に変換するには、v(t, X) = u(T - t, X)を設定して時間反転を適用します。

symsv(t, X)eq2 = subs(eq, {u, t}, {v(T - t, X), T - t}); eq2 = feval(symengine,'rewrite', eq2,'diff')

eq2 =
                  
                    $\frac{{σ (T - t, X)}^{2} \frac{\partial^{2}}{\partial X^{2}} v (t, X)}{2} + μ (T - t, X) \frac{\partial}{\partial X} v (t, X) - \frac{\partial}{\partial t} v (t, X) - v (t, X) r (T - t, X)$

ソルバーpdepeは次の形式の偏微分方程式を必要とします。ここで、係数c、fおよびsは、x、t、vおよび $\partial v / \partial x$ の関数です。

$c \frac{\partial v}{\partial t} = \frac{\partial f}{\partial x} + s$

方程式eq2をソルバーpdepeで解くには、eq2をこの形式にマッピングします。最初に、eq2の係数と項を抽出します。

symsdvtdvXXeq3 = subs(eq2, {diff(v, t), diff(v, X, X)}, {dvt, dvXX}); [k,terms] = coeffs(eq3, [dvt, dvXX])

k =
                  
                    $(\begin{array}{c} - 1 & \frac{{σ (T - t, X)}^{2}}{2} & μ (T - t, X) \frac{\partial}{\partial X} v (t, X) - v (t, X) r (T - t, X) \end{array})$

terms =
                   $(\begin{array}{c} dvt & dvXX & 1 \end{array})$

次に、eq2をk(1)*terms(1) + k(2)*terms(2) + k(3)*terms(3) = 0として書き換えます。時間微分の項を方程式の左辺に移動すると、次が得られます。

$\frac{\partial v}{\partial t} = k (2) \frac{\partial^{2} v}{\partial X^{2}} + k (3)$

右辺を個々に手動で部分積分します。結果は次のようになります。

$\frac{\partial}{\partial X} (k (2) \frac{\partial v}{\partial X}) + k (3) - \frac{\partial v}{\partial X} \frac{\partial k (2)}{\partial X}$

したがって、eq2をpdepeに適した形式で記述するには、次のパラメーターを使用します。

c = 1; f = k(2) * diff(v, X); s = k(3) - diff(v, X) * diff(k(2), X);

2.伊藤ルールの適用

資産価格は乗算過程に従います。つまり、価格の対数は SDE で表せますが、利益を表すのは価格の期待値のため、価格の期待値自体が問題となります。そのため、後者の SDE が必要です。

一般的に、確率過程Xは SDE として与えられ、伊藤ルールによって変換過程G(t, X)が次を満たすことが示されます。

$d G = (μ \frac{d G}{d X} + \frac{σ^{2}}{2} \frac{d^{2} G}{d X^{2}} + \frac{d G}{d t}) d t + \frac{d G}{d X} σ d B (t)$

価格の対数は 1 次元加算ブラウン運動で与えられると仮定します。つまり、muとsigmaは定数です。伊藤ルールを適用する関数を定義し、それを使用して加算ブラウン運動を幾何ブラウン運動に変換します。

ito = @(mu, sigma, G, t, X)...deal( mu * diff(G, X) + sigma^2/2 * diff(G, X, X) + diff(G, t),...diff(G, X) * sigma ); symsmu0sigma0[mu1, sigma1] = ito(mu0, sigma0, exp(X), t, X)

mu1 =
                  
                    $\frac{e^{X} {σ_{0}}^{2}}{2} + μ_{0} e^{X}$

sigma1 =
                   $σ_{0} e^{X}$

exp(X)をYで置き換えます。

symsYmu1 = subs(mu1, X, log(Y)); sigma1 = subs(sigma1, X, log(Y)); f = f(t, log(Y)); s = s(t, log(Y));

簡単にするために、利率を 0 と仮定します。これは、コルモゴロフの後退方程式としても知られる特別なケースです。

r0 = 0; c = subs(c, {mu, sigma}, {mu1, sigma1}); s = subs(s, {mu, sigma, r}, {mu1, sigma1, r0}); f = subs(f, {mu, sigma}, {mu1, sigma1});

3.ファインマン・カッツ方程式の求解

シンボリック式を MATLAB 関数ハンドルに変換する前に、diff(v(t, X), X)やv(t, X)などの関数呼び出しを変数に置き換えなければなりません。任意の有効な MATLAB 変数名を使用できます。

symsdvdxV; dvX = diff(v, X); c = subs(c, {dvX, v}, {dvdx, V}); f = subs(f, {dvX, v}, {dvdx, V}); s = subs(s, {dvX, v}, {dvdx, V});

平坦形状 (並進対称性) の場合、次の値を設定します。

m = 0;

また,数値をシンボリックパラメーターに代入します。

muvalue = 0; sigmavalue = 1; c0 = subs(c, {mu0, sigma0}, {muvalue, sigmavalue}); f0 = subs(f, {mu0, sigma0}, {muvalue, sigmavalue}); s0 = subs(s, {mu0, sigma0}, {muvalue, sigmavalue});

matlabFunctionを使用して関数ハンドルを作成します。係数c0、f0およびs0をpdepeで必要な形式、つまり、3 つの出力引数をもつ関数ハンドルの形式で渡します。

FeynmanKacPde = matlabFunction(c0, f0, s0,'Vars', [t, Y, V, dvdx])

FeynmanKacPde =function_handle with value:@(t,Y,V,dvdx)deal(1.0,0.0,0.0)

最終条件として、恒等写像を採用します。つまり、時間 $T$ における払い戻しは、資産価格そのものによって与えられます。金融派生商品を調べるために、この行を変更できます。

FeynmanKacIC = @(x) x;

PDE系の数値的求解は有限領域にのみ適用できます。したがって、境界条件を指定しなければなりません。資産価格が一定の範囲を上回るまたは下回るときに資産を売却するので、解vは境界ポイントxでx - v = 0を満たすと仮定します。別の境界条件を選択できます。たとえば、v = 0を使用してノックアウトオプションをモデル化できます。2 番目と 4 番目の出力の 0 は、境界条件が $\frac{d v}{d x}$ に依存しないことを示します。

FeynmanKacBC = @(xl, vl, xr, vr, t)...deal(xl - vl, 0, xr - vr, 0);

空間グリッドを選択します。これは、価格xの値の範囲です。左境界を 0 に設定します。これは、幾何ランダムウォークでは到達できません。右境界を 1 に設定します。資産価格がこの制限値に到達すると、資産は売却されます。

xgrid = linspace(0, 1, 101);

時間グリッドを選択します。最初に時間反転を適用しているため、これは最終時間Tまでの残り時間を表します。

tgrid = linspace(0, 1, 21);

方程式を解きます。

sol = pdepe(m,FeynmanKacPde,FeynmanKacIC,FeynmanKacBC,xgrid,tgrid);

解をプロットします。推定売価は、時間tにおける価格にほぼ線形に依存し、tにもわずかに依存します。

surf(xgrid, tgrid, sol) title('Expected selling price of asset'); xlabel('price'); ylabel('time'); zlabel('expected final price');

Figure contains an axes object. The axes object with title Expected selling price of asset contains an object of type surface.

時間tにおけるドリフト $μ 1$ の幾何ブラウン運動の状態は、初期値の $\exp (μ_{1} t)$ 倍の期待値を持つ対数正規分布された確率変数です。これは、制限値への到達のために売却されることのない資産の推定売価を表します。

Expected = transpose(exp(tgrid./2)) * xgrid;

上記で得られた解を期待値で除算することで、制限値に達する前に売却した場合に失われる利益が示されます。

sol2 = sol./Expected; surf(xgrid, tgrid, sol2) title('Ratio of expected final prices: with / without sales order at x=1') xlabel('price'); ylabel('time'); zlabel('ratio of final prices');

Figure contains an axes object. The axes object with title Ratio of expected final prices: with / without sales order at x=1 contains an object of type surface.

たとえば、売却注文に制限値が設定されている資産の期待支払い率と売却注文制限のない同じ資産の期待支払い率をtの関数として時間範囲T上にプロットします。注文の制限値が現在の価格の 2 倍のケースと 4 倍のケースをそれぞれ考えます。

plot(tgrid, sol2(:, xgrid == 1/2)); holdonplot(tgrid, sol2(:, xgrid == 1/4)); legend('Limit at two times current price','Limit at four times current price'); xlabel('timespan the order is valid'); ylabel('ratio of final prices: with / without sales order'); holdoff

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Limit at two times current price, Limit at four times current price.

4.期待される資産売却の時間の計算

教科書どおりの結果では、制限値に到達した場合と資産が売却された場合の期待終了時間は次の方程式で与えられます。

symsy(X)exitTimeEquation(X) = subs(eq, {r, u}, {0, y(X)}) == -1

exitTimeEquation(X) =
                  
                    $\frac{{σ (t, X)}^{2} \frac{\partial^{2}}{\partial X^{2}} y (X)}{2} + μ (t, X) \frac{\partial}{\partial X} y (X) = - 1$

これに加えて、yは境界において 0 でなければなりません。muとsigmaに対し、検討対象の実際の確率過程を挿入します。

exitTimeGBM = subs(subs(exitTimeEquation, {mu, sigma}, {mu1, sigma1}), Y, X)

exitTimeGBM(X) =
                  
                    $(\frac{X {σ_{0}}^{2}}{2} + X μ_{0}) \frac{\partial}{\partial X} y (X) + \frac{X^{2} {σ_{0}}^{2} \frac{\partial^{2}}{\partial X^{2}} y (X)}{2} = - 1$

この問題を任意の区間境界aおよびbについて解きます。

symsabexitTime = dsolve(exitTimeGBM, y(a) == 0, y(b) == 0)

exitTime =
                  
                    $\begin{array}{l} \frac{2 μ_{0} σ_{5} \log (b) - 2 μ_{0} σ_{4} \log (a) + a^{σ_{1}} {σ_{0}}^{2} σ_{5} σ_{4} - b^{σ_{1}} {σ_{0}}^{2} σ_{5} σ_{4}}{σ_{3}} - \frac{\log (X)}{μ_{0}} + \frac{σ_{2} (σ_{7} - σ_{6} - a^{σ_{1}} {σ_{0}}^{2} σ_{5} + b^{σ_{1}} {σ_{0}}^{2} σ_{4})}{σ_{3}} + \frac{X^{σ_{1}} {σ_{0}}^{2} σ_{2}}{2 {μ_{0}}^{2}} \\ where \\ σ_{1} = \frac{2 μ_{0}}{{σ_{0}}^{2}} \\ σ_{2} = e^{- \frac{2 μ_{0} \log (X)}{{σ_{0}}^{2}}} \\ σ_{3} = 2 {μ_{0}}^{2} (σ_{5} - σ_{4}) \\ σ_{4} = e^{- \frac{σ_{6}}{{σ_{0}}^{2}}} \\ σ_{5} = e^{- \frac{σ_{7}}{{σ_{0}}^{2}}} \\ σ_{6} = 2 μ_{0} \log (b) \\ σ_{7} = 2 μ_{0} \log (a) \end{array}$

mu0 = 0を直接代入することはできないので、代わりにmu0 -> 0の極限値を計算します。

L = limit(subs(exitTime, sigma0, sigmavalue), mu0, muvalue)

L =
                   $- (\log (X) - \log (a)) (\log (X) - \log (b))$

Lはa = 0では未定義です。0 < X < 1という仮定を設定します。

assume(0 < X & X < 1)

右境界で値b = 1を使用して、極限値を計算します。

limit(subs(L, b, 1), a, 0,'Right')

ans =
                   $\infty$

ファインマン・カッツの公式による確率過程のシミュレーション

1.確率微分方程式を使用したモデルのパラメーターの定義

2.伊藤ルールの適用

3.ファインマン・カッツ方程式の求解

4.期待される資産売却の時間の計算

期待される終了時間は無限大です。