ヨーロピアンオプションとインプライドボラティリティの数値計算

目的

ヨーロピアンオプションの価格をブラック・ショールズ式を用いて計算したり、ヨーロピアンオプションの価格からインプライドボラティリティを計算したりするプログラムを書いてみました。

デリバティブの数値計算をやるなら、すぐに出来ないと話にならないくらいの計算ですが、プログラミングに慣れていない人は、なかなか苦戦するようなので、ここでプログラム例を紹介しておきます。

正規分布の累積密度関数

ヨーロピアンオプションの価格を計算するには標準正規分布の累積密度関数が必要になります。コンパイラによっては誤差関数 erf() を用いて計算することも出来るようですが、対応していないコンパイラも多いようですので、自前で実装することにします。こちらで紹介した誤差関数を計算するアルゴリズムを書き換えて実現します。

#include<math.h>

double NormalCDF(double x){
    int n;
    double a, c, z;
    z = fabs(x);
    if(z<=3.5){
        c = 1;
        a = 0;
        for(n=1; n<=44; n++){
            a += x/(2*n-1)*c;
            c *= -x*x/(2*n);
        }
        return 0.5+a/sqrt(2*M_PI);
    } else {
        a = 0;
        if(z<1.0E50){
            for(n =44; n>=1; n--){
                a = n/(z+a);
            }
            a = exp(-z*z/2)/(z+a)/sqrt(2*M_PI);
        }
        if(x<0){
            return a;
        } else {
            return 1-a;
        }
    }
}

多倍長演算を用いて確かめたところ、どの x に対しても相対誤差は 8.0×10^-18 以下に抑えられています。これは、倍精度実数を用いる上では十分な精度です。

オプション価格の計算

標準正規分布の累積密度関数が用意できれば、ブラック・ショールズ式により、ヨーロピアンオプションの価格を計算することができます。

#include<math.h>

double NormalCDF(double x){
    中略
}

double EuropeanCall(double S, double K, double r, double v, double T, double t){
    double d1, d2, a, b;
    a = log(S/K);
    b = v*sqrt(T-t);
    d1 = (a+(r+v*v/2)*(T-t))/b;
    d2 = (a+(r-v*v/2)*(T-t))/b;
    return S*NormalCDF(d1) - K*exp(-r*(T-t))*NormalCDF(d2);
}

double EuropeanPut(double S, double K, double r, double v, double T, double t){
    double d1, d2, a, b;
    a = log(S/K);
    b = v*sqrt(T-t);
    d1 = (a+(r+v*v/2)*(T-t))/b;
    d2 = (a+(r-v*v/2)*(T-t))/b;
    return K*exp(-r*(T-t))*NormalCDF(-d2) - S*NormalCDF(-d1);
}

ここで、Sは原市場の現在価格、Kは権利行使価格、rは無リスク金利、vはボラティリティ、Tは満期、tは現在時間、です。

インプライドボラティリティの計算

オプション価格からボラティリティを逆算する方法について示します。ニュートン法を用いて非線形方程式を解くのが一番スマートかもしれませんが、ここでは応用範囲の広さを優先し、デリバティブの価格がボラティリティの増加関数である、ということだけを仮定します。この場合、微分可能性を仮定しないのでニュートン法は使えず、必然的に二分法を用いることになります。プログラムは以下のようになります。

#include<math.h>
#include<stdio.h>

double NormalCDF(double x){
    中略
}

double EuropeanCall(double S, double K, double r, double v, double T, double t){
    中略
}

double EuropeanPut(double S, double K, double r, double v, double T, double t){
    中略
}

#define Price(v) EuropeanCall(S,K,r,v,T,t)

double ImpliedVol(int *err, double P, double S, double K, double r, double v, double T, double t){
    double vh, vl, vm, w;

    for(vh = v; Price(vh)<P; vh*=2){
        if(vh > 1000){
            *err = 1;
            return 0.0;
        }
    }

    for(vl = v; Price(vl)>P; vl/=2){
        if(vl < 0.001){
            *err = -1;
            return 0.0;
        }
    }

    w = vh - vl;
    while(w>1.0E-16){
        vm = (vh+vl)/2;
        if(Price(vm) > P){
            vh = vm;
        } else {
            vl = vm;
        }
        w /= 2;
    }
    *err = 0;
    return (vh+vl)/2;
}

int main(){
    int err;
    double vol;
    
    vol = ImpliedVol(&err, 1500, 19300, 19000, 0.01, 0.3, 1.0, 0.7);
    if(err>0){
        printf("ボラティリティを最大限に大きくしても合致する値が見つかりません。\n",vol);
        return 1;
    }
    if(err<0){
        printf("ボラティリティを最小限に小さくしても合致する値が見つかりません。\n",vol);
        return 1;
    }
    printf("インプライドボラティリティは %.15f です。\n",vol);
    return 0;
}

枠をはみ出しているのは大目に見て下さい。「#define Price(v) …」の部分を書き換えれば、Putオプションの計算もできます。

main 関数から ImpliedVol(&err, 1500, 19300, 19000, 0.01, 0.3, 1.0, 0.7) を呼び出していますが、これは、オプション価格1500円、現在価格19300円、権利行使価格19000円、無リスク金利1％、期間1年、0.7年経過したオプションのインプライドボラティリティを求めています。ここで6番目の引数の0.3は基準ボラティリティーです。これは適当な値で構いません。errにはエラー情報が返ります。

ImpliedVol 関数の中では、まずは、基準ボラティリティーから始めてボラティリティーの値を倍々に増やしながら、目的価格よりオプション価格が大きくなるボラティリティーを探します。次に、基準ボラティリティーから始めてボラティリティーの値を半分々々に減らしながら、目的価格よりオプション価格が小さくなるボラティリティーを探します。ボラティリティの上限と下限が求まったならば、後は2分法で上限と下限の値を狭めていきます。上限と下限の差が 10^-16 以下になれば終了します。

基準ボラティリティーの1000倍までに上限が、0.001倍までに下限が見つからなければ、エラーとなります。このプログラムを他のオプション価格を求めるために流用した場合、プログラムによっては、ボラティリティとしてあまりに大き過ぎる値や小さ過ぎる値を入力すると不安定になることもあります。その場合は、上限や下限である 1000 や 0.001 の値を変更して下さい。