2

私は、Matlab (2012b) の Black-Scholes 式を使用してインプライド ボラティリティを計算しようとしていますが、何らかの形で行使価格に問題があります。たとえば、blsimpv(1558,1440,0.0024,(1/12),116.4) は NaN を返します。おそらく関数に問題があると思ったので、インターネットで他のmatlabスクリプトを検索し、カスタムニーズに合わせてカスタマイズしましたが、残念ながらまだ有効なインプライドボラティリティを返すことができません.

function sigma=impvol(C,S,K,r,T)

   %F=S*exp((r).*T);
   %G=C.*exp(r.*T);

   %alpha=log(F./K)./sqrt(T);
   %beta=0.5*sqrt(T);

   %a=beta.*(F+K);
   %b=sqrt(2*pi)*(0.5*(F-K)-G);
   %c=alpha.*(F-K);

   %disc=max(0,b.^2-4*a.*c);
   %sigma0=(-b+sqrt(disc))./(2*a);

   i=-1000;
   while i<=5000
      sigma0=i/1000;
      sigma=NewtonMethod(sigma0);
      if sigma<=10 && sigma>=-10
          fprintf('This is sigma %f',sigma)
      end
      i=i+1;
    end
end

function s1=NewtonMethod(s0)

    s1=s0;
    count=0;
    f=@(x) call(S,K,r,x,T)-C;
    fprime=@(x) call_vega(S,K,r,x,T);

    max_count=1e4;

    while max(abs(f(s1)))>1e-7 && count<max_count
        count=count+1;
        s0=s1;
        s1=s0-f(s0)./fprime(s0);

    end

end

end

function d=d1(S,K,r,sigma,T)
   d=(log(S./K)+(r+sigma.^2*0.5).*(T))./(sigma.*sqrt(T));
end

function d=d2(S,K,r,sigma,T)
   d=(log(S./K)+(r-sigma.^2*0.5).*(T))./(sigma.*sqrt(T));
end

function p=Phi(x)
  p=0.5*(1.+erf(x/sqrt(2)));
end

function p=PhiPrime(x)
   p=exp(-0.5*x.^2)/sqrt(2*pi);
end

function c=call(S,K,r,sigma,T)
  c=S.*Phi(d1(S,K,r,sigma,T))-K.*exp(-r.*(T)).*Phi(d2(S,K,r,sigma,T));
end

function v=call_vega(S,K,r,sigma,T)
   v=S.*PhiPrime(d1(S,K,r,sigma,T)).*sqrt(T);
end

ただし、impvol(116.4,1558,1440,0.0024,(1/12)) を実行すると、残念ながら値 'Inf' が返されます。ニュートン・ラプソン法が収束しないという問題がありますが、これを解決する方法がわかりません。この問題を解決する方法を知っている人、またはインプライド ボラティリティを計算する方法を知っている人はいますか?

あなたの助けを前もってありがとう！敬具、

ヘンク

4

3 に答える 3

1

ニュートン ラプソン法は、インプライド ボラティリティではうまく機能しません。二分法を使用する必要があります (Matlab でどのように使用されているかはわかりません)。http://en.wikipedia.org/wiki/Bisection_methodで説明されています。完全を期すために、次のように機能します。

1) high=200%/year のような任意の高い (不可能な) ボラティリティを選択します。

2) 可能な限り低いボラティリティ (low=0%) を選択します。

2a) ボラティリティが 0% の場合のオプション プレミアムを計算します。実際のプレミアムがそれよりも低い場合、それはボラティリティがマイナスであることを意味します (これは「不可能」です)。

3) インプライド ボラティリティが検出されない場合:

3.1) 「高値」と「安値」が非常に近い場合 (たとえば、小数点以下 5 桁まで等しい場合)、いずれかがインプライド ボラティリティです。そうでなければ...

3.2) 「高い」と「低い」の平均を計算します。平均=(高+低)/2

3.3) 平均ボラティリティのオプション プレミアムを計算します。

3.4) 実際のプレミアムが p(avg) よりも高い場合、インプライド ボラティリティは avg と max の間にある必要があるため、min=avg とします。

3.4a) 実際のプレミアムが p(avg) よりも低い場合、インプライド ボリュームは min と avg の間になければならないため、max=avg とします。

bisect の主な欠点は、最大値を選択する必要があるため、関数はそれよりも大きなインプライド ボラティリティを見つけられないことです。しかし、実際の使用には 200%/年程度で十分です。

ベガは導関数であるため、ニュートンの方法に似たさらに別の方法を使用します。したがって、範囲に限定されませんが、小さなベガによるハンティングと失敗を回避するために「線形化」修正を加えています。

def implied_volatility(type, premium, S, K, r, s_dummy, t):
    if S <= 0.000001 or K <= 0.000001 or t <= 0.000001 or premium <= 0.000001:
        return 0.0

    s = 0.35

    for cycle in range(0, 120):
        ext_premium = type(S, K, r, s, t)
        if abs(premium - ext_premium) < 0.005:
            return s
        ext_vega = type.vega(S, K, r, s, t)
        # print S, K, r, s, t, premium, ext_premium, ext_vega
        if ext_vega < 0.0000001:
            # Avoids zero division if stuck 
            ext_vega = 0.0000001

        s_new = s - (ext_premium - premium) / ext_vega
        if abs(s_new - s) > 1.00:
            # estimated s is too different from previous;
            # it is better to go linearly, since
            # vega is too small to give a good guess
            if s_new > s:
                s += 1.0
            else:
                s -= 1.0
        else:
            s = s_new

        if s < 0.0:
            # No volatility < 0%
            s = 0.0001
        if s > 99.99:
            # No point calculating volatilities > 9999%/year
            return 100.0

    return 0.0

それでも、バイセクトが最善の策だと思います。

于 2014-02-05T19:18:38.863 に答える
0

blsimpv の出力が NaN の場合に試行錯誤のような計算を行う単純な関数を作成しました。これにより、計算時間が大幅に遅くなりますが、常に望ましい結果が得られます。

関数は以下に使用されるように示されています

BSIVC(t,i)= blsimpv(S(t,i),K,r,tau(t),HestonCiter(t,i))
 if isnan(BSIVC(t,i));
  BSIVC(t,i)= secondIVcalc(HestonCiter(t,i),S(t,i),K,r,q,tau(t))
 end

関数自体は次のように説明されています。

function IV= secondIVcalc(HestonC,S,K,r,q,T)
lowCdif = 1;
a=0;
while lowCdif>0.0001
a= a+0.00001
lowCdif = HestonC - BSCprice(S,K,r,q,a,T);
end
IV= a;
end

BSCprice は matlab の組み込み関数ではないことに注意してください。

コードをわかりやすくするために、BSCprice の形式は BSCprice (原資産価格、権利行使価格、金利、配当利回り、インプライド ボリューム、満期までの時間) です。

于 2017-12-29T13:09:30.043 に答える