FinanceとHPCの備忘録

Quantitative finance/High Performance Computingに関して学んだことのメモです。Derivative Pricing/HPC/CUDA/OpenCL

CMSとコンベクシティ調整について

今回はCMSの評価によく使用されるReplication methodについて。
原論文は以下の通り。
https://www.deriscope.com/docs/Hagan_Convexity_Conundrums.pdf

CMSはスワップレートを原資産とするため、これをLIBORを原資産とするキャップやフロアと同様にフォワード測度で評価しようとすると、スワップ測度でマルチンゲールになるスワップレートにドリフト項が生じる。そのため、CMSの評価にはこのドリフト項を調整するコンベクシティ調整を行う必要がある。
このコンベクシティ調整には以下の2通りの手法が存在する。

  • ATMボラを使用した解析解
  • スワップションによる静的Replication手法

ATMボラによる調整は計算が非常に軽くジョンハルのフィナンシャルエンジニアリングにも掲載されている古い手法である。対してReplicationによる調整は多くのスワップションを評価する必要があるため計算が重いが、ボラティリティスマイルを勘案することが可能で実務上でメジャーな手法となっている。そのため本記事ではReplication手法をメインに紹介する。

スワップ取引の無裁定条件より、スワップレートR_s(t)は割引債DFを使用して以下のように書ける。
R_s(t)=\frac{DF(t;s_0)-DF(t;s_n)}{A(t)}
L(t)=\sum_{j=1}^n \alpha_(j)DF(t;s_j)
ここで、s_0,s_nはスワップの開始日、最終日を示し、A(t)はアニュイティと呼ばれる。
スワップの満期が\tauでCMSキャップレットのペイオフ支払時期がt_pのとき、ペイオフは(R_s(\tau)-K))^{+}DF(\tau;t_p)となる。
よって求めたいPVは
V(0)=L(0)E[(R_s(\tau)-K))^{+}DF(\tau;t_p)]
ここではフィルトレーションは省略する。
アニュイティがニュメレールのときにマルチンゲールとなるようなスワップ測度のもとでは下記が満たされる。
E [DF(\tau;t_p)/L(\tau)]=DF(0;t_p)/L(0)
よって
V(0)=DF(t_p)E[(R_s(\tau)-K))^{+}\frac{DF(\tau;t_p)/L(\tau)}{DF(0;t_p)/L(0)}]
より、
V(0)=DF(t_p)E[(R_s(\tau)-K))^{+}]+DF(t_p)E[max(R_s(\tau)-K,0)(\frac{DF(\tau;t_p)/L(\tau)}{DF(0;t_p)/L(0)}-1)]
となる。
ここで、第１項は解析的に解くことができ、第２項はコンベクシティ調整項となる。

続いてコンベクシティ調整項の評価のためには、DFとアニュイティの比を何らかパラメトリックな形に表現する必要がある。ここでEURでは一般的なCashで決済が行われるCash-settled-Swaptionを考えると、この形式では本来ゼロレートから計算される割引率がスワップレートによって近似的に計算されるものとして定義される。よって
DF(\tau;t_p)/L(\tau)=G(R_s(\tau))
G(R_s)=\frac{R_s}{(1+R_s/q)^\delta}\frac{1}{1-\frac{1}{(1+R_s/q)^n}}
となる。

次に静的Replicationの公式より
f'(K)[R_s-K]^{+}+\int_{K}^{\infty}[R_s-x]^{+}f''(x)dx=f(R_s)1_{R_s>K}
f(x)=[x-K](\frac{G(x)}{G(R_s^0)}-1)
これにより、CMSのコンベクシティ調整項をスワップションのストライクでの積分で表現することができた。

最後にG(x)を１次のテイラー展開することで、f(x)の2階微分が一定値となるので、Caplet全体の価値は
V(0)=DF(0;t_p)/L_0 C(K)+G'(R_s^0)[(K-R_s^0)C(K)+2\int_{K}^{\infty}C(x)dx]
ここでC(K)はスワップションの価値である。

積分項により、ストライク全体のスワップションを評価するため、ボラティリティのスキューやスマイルをSVやLVを使用して反映させることが可能となる。

まとめ
このReplication手法によりスキュー・スマイルを勘案したコンベクシティ調整が可能となり、実務では広く取り入れられた手法である。
ただし、実装においては下記の課題があり、様々な対処が行われている。
・高ストライクにおけるスワップション評価の安定性
・数値積分の計算コスト

論文 - 『A STOCHASTIC VOLATILITY FORWARD LIBOR MODEL WITH A TERM STRUCTURE OF VOLATILITY SMILES』

V.Piterbergによる2003年の論文を読んでみた。

題名の通り、フォワードLIBORにボラティリティの確率変動を勘案したモデル(FL-SV model)が紹介されている。特徴はTime DependentなVolatility SmileのSkewパラメータを導入することでSkewのダイナミクスを表現している点と、CalibrationにParameter Averagingを用いることだと思う。

例えば、CallOptionのPriceを求めるのに複雑な数値計算を要する場合はCalibration targetの数だけ繰り返して数値計算をしなければならない。しかし、Parameter Averagingを用いれば１回の計算のみで済むことから、計算コストに優れた手法となっている。

この論文は2003年とやや古いがこのParameter Averaging手法は様々なSV モデルに応用されており今でも読む価値は高いっぽい。以下内容のメモ。

 

dS(t)=\sigma(bS(t)+(1-b)S(0))\sqrt{z(t)}dW(t)

dz(t)=\theta(z_0-z(t))dt+\gamma\sqrt{z(t)}dV(t)

z(0)=z_0

dVdW=0

bはVolatility Skew、\thetaは平均回帰スピードを表現している。ただしS(t)はSwap Rate。Local VolatilityをbS(t)+(1-b)S(0)とすることで、b=0のときは対数正規分布モデル、b=0のときは正規分布モデルとなり、調整することで両者の中間分布を表現することができる。また、X_t=bS(t)+(1-b)S(0)とすればX_tを原資産とするCall Optionの解析解を求めることも可能なため、解析解がやや扱いにくいCEVに比してよく使われている。

Brownian Motionの相関が0なのが微妙だが、下記論文では\rho\not=0でかつVolatility SmileのConvexityパラメータも取り入れて拡張を行っている。

papers.ssrn.com

 

 Forward LIBORを割引債P(t,T)で表すと,

L_n(t)=\displaystyle \frac{P(_t,T_n)-P(t,T_{n+1})}{\tau_{n+1}P(t,T_{n+1})}

 

Swap RateはForward LIBORの集合で表現できるのでForward LIBORのSDEの確率変動を全部集めれば、無裁定条件下でSwap Rateの動きも全部説明ができる。ここらへんは元のLMMと一緒。

 

続いてTime DependentなパラメータをParameter Averagingを使って推定する。それにはまず以下のようなTime Dependentパラメータを持つSDEを新たに定義する。以下の式が与えられた下で、"Effective parameter"（Time Dependentパラメータの平均みたいなイメージ）はどうなるかを調べたいとする。

 

dS(t)=\sigma(t)(\beta(t)S(t)+(1-\beta(t))S(0))\sqrt{z(t)}dW(t)

 

つまり、Time Dependentパラメータが分かっていれば、同じCall Option価格となるようなConstantパラメータも分かるので、逆にConstantパラメータが分かればTImeDependentパラメータも分かる。

Calibration TargetはCall Optionを仮定する。冒頭で述べたように本来Time Dependentパラメータを推定する際は、Time Dependentモデルから得たCall Optionの価格を、Marketから得られたCall Optionの価格に最適化するため、Call Optionのモデル価格を数値計算等を用いて解かなければならない(Time Dependentの場合解析解が得られることはないため)。

しかし、Parameter Averagingを用いる場合、まず第一に

  • Constant Parameterモデルから得たCall Optionの価格をMarketから得られたCall Optionの価格にCalibrationする。その際、モデル価格は解析的に解けるため(Black-ScholesやHestonなど)、計算コストが低くなる
  • 求められたConstant Parameterから、同じCall Option価格になるようなTime Dependent Parameterを推定する。Time Dependent Parameterを持つモデルのCall Option価格を求めることなく、Calibrationを行うことができる。

つまり、CalibrationにConstant Parameterモデルを経由するかしないかの違いであるが、計算コストが格段に改善するためよく使用される。

 

次に気になるのが、Constant ParameterからTime Dependent Parameterを推定する方法。まず、ボラティリティ\sigmaを推定する。

  • ATMCall Option価格のBlack Scholesの公式をg(x)\approx a+b\exp(-cx)と近似することで扱いを楽にする
  • Time Dpendentな場合解析解が存在しないので、Call Option価格をラプラス変換することで特性関数を導く
  • ２つのCall Option価格が等しくなるように式変形を行うと、以下のようなTime DependentなボラティリティとConstantなボラティリティの関係式が導かれるので解析解・数値計算を使用してCalibrationを行う

E exp(-c\int_{0}^{T}\sigma^2z(t)dt) =E exp(-c\overline{\sigma}^2\int_{0}^{T}z(t)dt)

続いてSkewパラメータbを推定する。それには2つのSDEをSmall Noise Expansionしてその差分をとることで、\mathcal{O}(\epsilon^2)かそれ以下に抑えることができる。(Small Noise Expansionは誤解を恐れずに言えば例えば、Volatilityパラメータが\epsilon V(t)みたいに表現できて\epsilonがめちゃくちゃ小さいときは\epsilonについてTaylor ExpansionすることでUndelyingのSDE解を漸近的に求めることができる手法）

 

実装する際にはTime Dependent ParameterはPiecewise Constantなものとして扱う。

 

 

【CUDA】Monte Carlo Simulationを用いたOption Pricing

別記事で「Cooperative Group」を用いたReductionについて記事を書いたけど、肝心のMonteCarloについては触れていなかったので簡単にメモ。
gomiwo.hatenablog.com

モンテカルロでのプライシング理論はググれば沢山情報が出てくるのでそちらを参照。
簡単に書くと、UndelyingのSDE使って満期までのパスを生成→パスごとにペイオフを算出→パス本数で割る、という流れ。

以下はCUDAのサンプルに入っているMonteCarlo Simulationのサンプルコードを抜粋。

・求めたいOprion Priceが複数ある状況で、１つのOptionにつき１Block割り当て、１つの株価パスにつき1Thread割り当てている。
・「cooperative_group.h」を使ってBlock内のThreadを分割。Reductionで使用。
・Call Option Priceの配列はShared Memoryに格納
・株価パスの生成は#pragma unroll 8でループ展開することで速度向上

static __global__ void MonteCarloOneBlockPerOption(
    curandState * __restrict rngStates,
    const __TOptionData * __restrict d_OptionData,
    __TOptionValue * __restrict d_CallValue,
    int pathN,
    int optionN)
{
    // Handle to thread block group
    cg::thread_block cta = cg::this_thread_block();
    cg::thread_block_tile<32> tile32 = cg::tiled_partition<32>(cta); //Thread Groupを作成

    const int SUM_N = THREAD_N;
    __shared__ real s_SumCall[SUM_N];
    __shared__ real s_Sum2Call[SUM_N];

    // determine global thread id
    int tid = threadIdx.x + blockIdx.x * blockDim.x;

    // Copy random number state to local memory for efficiency
    curandState localState = rngStates[tid];
    for(int optionIndex = blockIdx.x; optionIndex < optionN; optionIndex += gridDim.x) //Block1つにつきOption Price1つを計算
    {
        const real        S = d_OptionData[optionIndex].S;
        const real        X = d_OptionData[optionIndex].X;
        const real    MuByT = d_OptionData[optionIndex].MuByT;
        const real VBySqrtT = d_OptionData[optionIndex].VBySqrtT;

        //Cycle through the entire samples array:
        //derive end stock price for each path
        //accumulate partial integrals into intermediate shared memory buffer
        for (int iSum = threadIdx.x; iSum < SUM_N; iSum += blockDim.x) //Thread１つにつきStockのPath1本を生成
        {
            __TOptionValue sumCall = {0, 0};

            #pragma unroll 8 //ループ展開で速度向上
            for (int i = iSum; i < pathN; i += SUM_N)
            {
                real              r = curand_normal(&localState); //正規乱数の生成
                real      callValue = endCallValue(S, X, r, MuByT, VBySqrtT); //PayOffの計算
                sumCall.Expected   += callValue;
                sumCall.Confidence += callValue * callValue;
            }

            s_SumCall[iSum]  = sumCall.Expected;
            s_Sum2Call[iSum] = sumCall.Confidence;
        }

        //Reduce shared memory accumulators
        //and write final result to global memory
        cg::sync(cta);
        sumReduce<real, SUM_N, THREAD_N>(s_SumCall, s_Sum2Call, cta, tile32, &d_CallValue[optionIndex]); //Reductionの実行
    }
}

__device__ inline float endCallValue(float S, float X, float r, float MuByT, float VBySqrtT)
{
    float callValue = S * __expf(MuByT + VBySqrtT * r) - X;
    return (callValue > 0.0F) ? callValue : 0.0F;
}

CUDAでMonte Carlo Simulation -「Cooperative Groups」

Cooperative Groupsについて学んだのでメモ。

詳細は以下参照。

Cooperative Groups: Flexible CUDA Thread Programming | NVIDIA Developer Blog
Cooperative Groups: Flexible CUDA Thread Programming | NVIDIA Developer Blog

 
普通にCUDAでスレッド並列しようとするとブロック単位で動かすことになってしまうけど、cooperative_group.hを使うことで、より粒度を細かくすることができるという代物。
BlockとThreadの間に一つ階層が増えるイメージだと思ってる。

以下はCUDAのサンプルにあるMonteCarlo SimulationのReductionのサンプルコード。
サンプルコードでは２つの配列を別々にReductionしてるが、１つのOption Priceを求めたいときは１つの配列で十分。

・this_thread_block()でBlock内のスレッドを扱うインスタンスを召喚
・tiled_partition<32>()でスレッドを32個ごとに分割する
・thread_rank()でthreadのidを指定
・sync()でスレッドを同期

・下記コードで用いられているReduction手法はインタリーブペア方式と呼ばれるもので、配列を半分に分けて分割した１要素目同士、２要素目同士…を足し合わせていく。２ループ目では配列をさらに半分に分割して同様に足し合わせていく。そして最終的に合計値が求まる。

#include <cooperative_groups.h>

namespace cg = cooperative_groups;

template<class T, int SUM_N, int blockSize>
__device__ void sumReduce(T *sum, T *sum2, cg::thread_block &cta, cg::thread_block_tile<32> &tile32, __TOptionValue *d_CallValue)
{
const int VEC = 32; //Threadのサイズ
const int tid = cta.thread_rank();

T beta = sum[tid];
T beta2 = sum2[tid];
T temp, temp2;

for (int i = VEC/2; i > 0; i>>=1)  //sum[tid]とsum[tid+i]を足し合わせる。１ループ目でVEC/2回足し合わせ、２ループ目でVEC/2/2回足し合わせていく
{
if (tile32.thread_rank() < i)
{
temp = sum[tid+i];
temp2 = sum2[tid+i];
beta += temp;
beta2 += temp2;
sum[tid] = beta;
sum2[tid] = beta2;
}
cg::sync(tile32); //分割したPartition Threadを同期
}
cg::sync(cta); //Block内のThreadを同期

if (tid == 0) //最終結果
{
beta = 0;
beta2 = 0;
for (int i = 0; i < blockDim.x; i += VEC)  //Partition間で合計値を計算
{
beta += sum[i];
beta2 += sum2[i];
}
__TOptionValue t = {beta, beta2};
*d_CallValue = t;
}
cg::sync(cta);
}