ファインマン・カッツの公式を使用した確率過程のシミュレーション

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

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

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

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) dt)$ として計算されます。
$μ (t, X)$ と $σ (t, X)$ は、確率過程 X(t) のドリフトと拡散です。

syms mu(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 は初期条件を処理します。最終条件を初期条件に変換するには、u(t, X) = v(T - t, X) を設定して時間反転を適用します。

syms v(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 の係数と項を抽出します。

syms dvt dvXX
eq3 = 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 );

syms mu0 sigma0
[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 で置き換えます。

syms Y
mu1 = 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 変数名を使用できます。

syms dvdx V;
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,(Y.^2.*dvdx)./2.0,(Y.*dvdx)./2.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, xlabel price, ylabel time 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, xlabel price, ylabel time contains an object of type surface.

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

plot(tgrid, sol2(:, xgrid == 1/2));
hold on
plot(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');
hold off

Figure contains an axes object. The axes object with xlabel timespan the order is valid, ylabel ratio of final prices: with / without sales order contains 2 objects of type line. These objects represent Limit at two times current price, Limit at four times current price.

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

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

syms y(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 について解きます。

syms a b
exitTime = dsolve(exitTimeGBM, y(a) == 0, y(b) == 0)

exitTime = 
 $\begin{array}{l} \frac{2 μ_{0} σ_{4} \log (b) - 2 μ_{0} σ_{5} \log (a) + a^{σ_{1}} {σ_{0}}^{2} σ_{4} σ_{5} - b^{σ_{1}} {σ_{0}}^{2} σ_{4} σ_{5}}{σ_{3}} - \frac{\log (X)}{μ_{0}} + \frac{σ_{2} (σ_{6} - σ_{7} - a^{σ_{1}} {σ_{0}}^{2} σ_{4} + b^{σ_{1}} {σ_{0}}^{2} σ_{5})}{σ_{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} (σ_{4} - σ_{5}) \\ σ_{4} = e^{- \frac{σ_{6}}{{σ_{0}}^{2}}} \\ σ_{5} = e^{- \frac{σ_{7}}{{σ_{0}}^{2}}} \\ σ_{6} = 2 μ_{0} \log (a) \\ σ_{7} = 2 μ_{0} \log (b) \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$

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

参考

関数

diff | subs | coeffs | matlabFunction | dsolve