G++ calibration （工事中）

〜
class G_n_plusplus:
    def __init__(self, n_factor = 2) -> None:
        """
        Augs
            discount_factor_vec: 割引率のベクトル。スワップのCFに合う形を想定
        """
        # ファクターの数を設定する。
        self.n_factor = n_factor        

        # 初期値を設定する
        self.variance_vec = np.array([1] * n_factor)
        self.drift_vec = np.array([1] * n_factor)
        self.a_vec = np.array([1] * n_factor)
        self.correlation_mat = np.eye(self.n_factor)


    def _calc_swaption_price(self, maturity, strike, nominal_value, first_swap_time, time_fraction, n_swap, discount_factor_vec) -> float:
        """
        Augs
            marurity: スワップション満期（年）
            strike: スワップションの行使価格
            nominal_value: よくわからんけど、論文だと"N"
            first_swap_time: スワップ行使の開始時点（年）
            time_fraction: スワップ行使の時間幅(年)
            n_swap: スワップ回数
            discount_factor_vec: 割引率のベクトル。スワップのCFに合う形を想定(長さは"n_swap"と一致。)
        """ 
        # 入力を整理
        self.maturity = maturity
        self.strike = strike
        self.nominal_value = nominal_value
        self.first_swap_time = first_swap_time
        self.time_fraction = time_fraction
        self.n_swap = n_swap
        self.discount_factor_vec = discount_factor_vec

        # 最後のスワップの時刻をfinal_swap_timeとする
        self.final_swap_time = self.first_swap_time + self.time_fraction * self.n_swap

        # time fractionとDiscount Factorの積和を計算する。（ここでは、time fractionが定数と仮定しているため、Discount Factorの和にtime fractionをあとから掛けている）
        self.P = sum(self.discount_factor_vec) * self.time_fraction
        
        # Aのリストを作成する
        self._calc_A()

        # VOLを計算する
        self._calc_VOL()
        
        # 理論解を計算
        self.theoretical_swaption_price = self.nominal_value * self.VOL * self.P / math.sqrt(2 * np.pi)

    def _calc_A(self) -> list:
        """
        A_i (i = 1, 2, ..., factor数)を作成する。
        calc_swaption_priceのサブ的な関数
        """
        self.A = []
        for i in range(self.n_factor):
            # A の二つ目の項と最後の項を分けて計算
            # self.discount_factor_vec[0]は"P(0, T)"を表す
            A = math.exp(- self.drift_vec[i] * self.maturity) * self.discount_factor_vec[0] / self.P - math.exp(- self.drift_vec[i] * self.final_swap_time) * self.discount_factor_vec[-1] / self.P

            final_term_of_A = 0
            for j in range(self.n_swap):
                # スワップの実行時刻
                t_j = self.first_swap_time + self.time_fraction * j

                final_term_of_A += math.exp(- self.drift_vec[i] * t_j) * self.time_fraction * self.discount_factor_vec[j]
            
            final_term_of_A *= - self.strike
            
            # 最後の項を加算
            A += final_term_of_A
            
            # リストに格納
            self.A.append(A)
        
    def _calc_VOL(self) -> float:
        """
        VOLを計算する。calc_swaption_priceのサブ的な関数
        """
        VOL = 0
        for i in range(self.n_factor):
            for j in range(self.n_factor):
                VOL += self.variance_vec[i] * self.variance_vec[j] * self.correlation_mat[i, j] * self.A[i] * self.A[j] * \
                    (math.exp((self.drift_vec[i] + self.drift_vec[j]) * self.maturity) - 1) / (self.drift_vec[i] + self.drift_vec[j])
                
        self.VOL = math.sqrt(VOL)

maturity = 10
strike = 0.05
nominal_value = 1
first_swap_time = 11
time_fraction = 0.5
n_swap = 6
discount_factor_vec = [0.99, 0.95, 0.85, 0.81, 0.8, 0.75]

plt.plot(discount_factor_vec)

instance_G_2 = G_n_plusplus()

instance_G_2._calc_swaption_price(maturity=maturity, strike=strike, nominal_value=nominal_value, first_swap_time=first_swap_time, time_fraction=time_fraction, n_swap=n_swap, discount_factor_vec=discount_factor_vec)
instance_G_2.theoretical_swaption_price

〜

キュムラント展開

メモです。

ある確率変数の積率母関数とキュムラント母関数は以下の様に定義される。

\begin{equation}
M_X(t) = E[e^{t X}]
\end{equation}

\begin{equation}
C_X(t) = \log M_X(t)
\end{equation}

従って、マクローリン展開を用いると以下の式が得られる。

\begin{eqnarray}
C_X(t) &= \sum^{\infty}_{n = 0} \frac{1}{n!} \frac{d^n C_X(t)}{d t^n} t^n \\
&= \sum^{\infty}_{n = 0} \frac{c_n}{n!} t^n \tag{a}
\end{eqnarray}
ここで、であり、これを次のキュムラントと呼ぶ。

この時、以下の式が成り立つ。

\begin{equation}
\label{eq:goal}
\begin{split}
・\nu_1 &= c_1\\
・\nu_2 &= c_1^2 + c_2\\
・\nu_3 &= c_1^3 + 3c_1 c_2 + c_3^3
\end{split}
\end{equation}
ここで、] であり、つまり次のモーメントである。

証明

定義より

\begin{equation}
\label{eq1}
\begin{split}
M(t) &= \int^{\infty}_{-\infty} e^{tx} f_X{t}dx \\
&= \int^{\infty}_{-\infty} \sum^{\infty}_{m = 0} \frac{(t x)^m}{m!}f_X(x) dx \\
&= \sum^{\infty}_{m = 0} \frac{t^m}{m!} \int^{\infty}_{-\infty} x^m f_X(x) dx \\
&= \sum^{\infty}_{m = 0} \frac{t^m}{m!} \nu_m
\end{split}
\end{equation}
（積分と無限級数の交換を行っている（若干怪しい））

また、より

\begin{eqnarray}
M_X(t) &= \sum^{\infty}_{m = 0} \frac{C_X(t)}{n!} \tag{b}
\end{eqnarray}
(a)式を代入する。また、であるので

\begin{eqnarray}
M_X(t) &= \sum^{\infty}_{m = 0} (\sum^{\infty}_{n = 1} \frac{c_n}{n!} t^n) \\
&= 1 + c_1 t + (\frac{1}{2} c_1^2 + \frac{1}{2} c_2) t^2 + (\frac{1}{6} c_1 ^ 3 + \frac{1}{2} ・2 c_1 ・\frac{1}{2} c_2 + \frac{1}{6}c_3) t^3 + \dots \\
&= 1 + c_1 t + \frac{1}{2!} (c_1^2 + c_2) t^2 + \frac{1}{3!}(c_1 ^ 3 + 3 c_1 c_2 + c_3) t^3 + \dots \tag{c}
\end{eqnarray}
となる。（無限級数の順番を変えてる、また怪しい）

従って、(b), (c)式で係数比較を行うと題意は示せる。

結果は間違いないのですが、証明が怪しさ満点になってしまいました。
時間があれば今度ちゃんと示そう。。

ブラウン運動の二乗の期待値に関する定理

いつか使うかもしれない（使わないかもしれない）定理に関するメモを載せる。

定理1：

を正の実数とし、をブラウン運動とする。この時、任意の自然数に対して
\begin{equation}
E[W_T^{2n}] = \frac{(2n)!}{2^n n!} T^n
\end{equation}
が成り立つ。
定理１を証明する前に、以下の定理２と補題を用意する。

定理２：
梶原壌二（2005）.「改訂増補　独修微分積分学」　現代数学社

のp.86から拝借する。

関数がで連続で（条件０）

\begin{equation}
\begin{split}
g(t) &= \int^{\infty}_0 f(x, t) dx\\
&= \lim_{n \rightarrow \infty} \int^{n}_0 f(x, t) dx
\end{split}
\end{equation}
は収束する（有限な収束先が存在）（条件１）。さらに、偏導関数が存在して連続であり（条件２）、と無関係な関数が存在して

\begin{equation}
\int^{\infty}_0 h(x) dx < \infty
\end{equation}
（条件３）と

\begin{equation}
|f_t(x, t)| \leq h(x)
\end{equation}
（条件４）が成り立つ時、

積分と微分の交換
\begin{equation}
g'(t) = \frac{d}{dt} \int^{\infty}_0 f(x, t) = \int^{\infty}_0 \frac{\partial}{\partial t} f(x,t) dx
\end{equation}
とすることが出来る。

（証明は長くなりそうなので逃げます！！（謝罪）
尚、類似の定理が以下のページに証明付きで載っています。以下の系で扱う関数については、こちらの定理の条件も問題なく満たしてる（はず））解析学基礎/関数列の極限 - Wikibooks

系：
に対して、
\begin{eqnarray}
g_k (a) &= \int^{\infty}_0 f_k(a,x) dx \\
&= \int^{\infty}_0 x^{2k} e^{-ax^2} dx \tag{a}
\end{eqnarray}
とする。この時、

\begin{eqnarray}
\frac{d}{da} g_k(a, x) &= \int^{\infty}_0 \frac{\partial}{\partial a} f_k (a, x) dx \\
&= - \int^\infty_0 x^{2(k + 1)} e^{-ax^2} dx \tag{b}
\end{eqnarray}
が成り立つ。

系の証明：
定理２が使える事を確認する。まず、はそれぞれ連続関数の和で表されているので連続である（条件０、２）。また()に関して、と置換を行うとであることから

\begin{equation}
\begin{split}
g_k (a) &= \int^{\infty}_0 (\sqrt{\frac{t}{a}})^{2k} e^{-t} (\frac{1}{2\sqrt{at}} dt) \\
&= \frac{1}{2 a^{k + \frac{1}{2}}} \int^{\infty}_0 t^{(k + \frac{1}{2}) - 1} e^{-t} dt \\
&= \frac{1}{2 a^{k + \frac{1}{2}}} \Gamma(k + \frac{1}{2})
\end{split}
\end{equation}
と計算する事が出来る。ただし、はガンマ関数である。従って、条件１も満たされる。更に、とすれば

\begin{equation}
\begin{split}
|\frac{\partial}{\partial a} f_k(a, x)| &= x^{2 (k + 1)} e^{- a x^2} \\
& \leq h_k (x)
\end{split}
\end{equation}
である（条件４）。ここで、について考える。と置換するとであるから、

\begin{equation}
\begin{split}
\int^{\infty}_0 h_k(x) dx &= \int^{\infty}_0 x^{2(k + 1)} e^{-x^2} dx \\
&= \int^{\infty}_0 t^{k + 1} e^{-t} (\frac{dt}{2 \sqrt{t}}) \\
&= \frac{1}{2} \int^{\infty}_0 t^{(k + \frac{1}{2}) - 1} e^{-t} dt \\
&= \frac{1}{2} \Gamma(k + \frac{1}{2}) < \infty
\end{split}
\end{equation}
が示せる（条件３）。

従って、定理２を使うための条件を示すことが出来た。(a)式に定理２を適用すれば、(b)式が得られる。

補題：

\begin{equation}
\label{eq:3-3}
\begin{split}
\int_ {-\infty} ^{\infty} x^{2n} e^{-ax^2} dx = (-1)^n \frac{d^n}{da^n} \left( \sqrt{\frac{\pi}{a}} \right) \\
\end{split}
\end{equation}

補題の証明：
これは、ガウス積分

\begin{eqnarray}
\int^{\infty}_{-\infty} e^{-ax^2} dx = \sqrt{\frac{\pi}{a}} \tag{c}
\end{eqnarray}
の両辺をで回微分すれば得られる。左辺について考えると、被積分関数は偶関数なので

\begin{equation}
\frac{d^n}{da^n} \int^{\infty}_{-\infty} e^{-ax^2} dx = 2 \frac{d^n}{da^n} \int^{\infty}_{0} e^{-ax^2} dx
\end{equation}
となり、これに上で示した系を適用することで、

\begin{equation}
\begin{split}
\frac{d^n}{da^n} \int^{\infty}_{-\infty} e^{-ax^2} dx &= 2 \int^{\infty}_{0} \frac{\partial^n}{\partial a^n} e^{-ax^2} dx \\
&= 2 (-1)^n \int_{0} ^{\infty} x^{2n} e^{-ax^2} dx \\
&= (-1)^n \int_ {-\infty} ^{\infty} x^{2n} e^{-ax^2} dx \\
\end{split}
\end{equation}
となり、(c)式は示せた。

定理１の証明：
なので、期待値の定義より
\begin{equation}
\label{eq:3-1}
\begin{split}
E[W_T^{2n}] = \int_ {-\infty} ^{\infty}x^{2n} \frac{1}{\sqrt{2 \pi T}} e^{-\frac{x^2}{2T}} dx
\end{split}
\end{equation}
という変数変換を用いて、

\begin{equation}
\label{eq:3-2}
\begin{split}
E[W_T^{2n}] &= (\sqrt{T} y)^{2n} \int_ {-\infty} ^{\infty} \frac{1}{\sqrt{2 \pi T}} e^{-\frac{y^2}{2}} (\sqrt{T} dy) \\
&= T^n \frac{1}{\sqrt{2 \pi}} \int_ {-\infty} ^{\infty} y^{2n} e^{-\frac{y^2}{2}} dy
\end{split}
\end{equation}
ここで補題を用いて、

\begin{equation}
\label{eq:3-4}
\begin{split}
E[W_T^{2n}] &= T^n \frac{1}{\sqrt{2 \pi}} \{ (-1)^n (\frac{d^n}{da^n} a^{-\frac{1}{2}} )_{|a =\frac{1}{2}} \} \\
&= \frac{T^n}{\sqrt{2}} \{ (-1)^n (\frac{d^n}{da^n} a^{-\frac{1}{2}} )_{|a =\frac{1}{2}} \}
\end{split}
\end{equation}
が得られる。さらにこれを計算すると、

\begin{equation}
\label{eq:3-5}
\begin{split}
E[W_T^{2n}] &= \frac{T^n} {\sqrt{2}} (-1)^n \{ \left( -\frac{1}{2} \right) \left( -\frac{3}{2} \right) ・・・　\left( -\frac{2n - 1}{2} \right) a ^{-\frac{1}{2} - n }　\}_{|a =\frac{1}{2}} \\
&= \frac{T^n} {\sqrt{2}} (-1)^n (-1)^n \{ \left( \frac{1}{2} \right) \left( \frac{3}{2} \right) ・・・　\left( \frac{2n - 1}{2} \right) 2 ^{\frac{1}{2} + n } \} \\
&= \frac{T^n} {\sqrt{2}} \frac{1・3・・・(2n-3)(2n-1)}{2^n} 2^{\frac{1}{2} + n} \\
&= T^n \{1・3・・・(2n-3)(2n-1) \}
\end{split}
\end{equation}
が得られる。これを階乗を用いて表現すると、

\begin{equation}
\label{eq:3-6}
\begin{split}
E[W_T^{2n}] = \frac{(2n)!}{2^n (n!)} T^n
\end{split}
\end{equation}
という結果が得られる。（証明終わり）

数学的に怪しい所とか、誤植等へのツッコミやご指摘は大歓迎です。コメント欄やTwitterを通じて連絡していただけると大変助かります。

参考文献

・
blog.goo.ne.jp
・
独修微分積分学 | 梶原 壌二 |本 | 通販 | Amazon

・
ja.wikibooks.org

オプショングリークスの計算

ブラックショールズモデルにおけるオプショングリークスの計算式とプログラムを載せるだけ
何故か日本語で詳しく計算を追ってる記事が見つからなかったため

グリークスとは

株式に配当が含まれている場合の、ブラックショールズ方程式におけるヨーロピアンコール/プットオプション価格()は現在時刻t、原資産価格、ボラティリティ、行使価格K、無リスク金利r、配当利回りq、満期Tまでの時間によって決まる。具体的には以下の式

ただし、

グリークスとはこの式を各変数で偏微分したものであり、つまりは各変数に対するオプション価格の感応度を表したものである。 具体的には以下の通り¹

計算

準備

この先度々出てくる式を先に計算する。

１.（微積分の交換）

２.

３.

４.（プットコールパリティ）

デルタ（Δ）

コール/プットオプションのデルタをそれぞれと表す。

コールのデルタ

ここで、

と(1)式より

が得られる。更に二項目は(2)式を用いることで、

であることから０である。従って、

プットのデルタ

(4)式より、

よって、

ガンマ（γ）

コール/プットオプションのデルタをそれぞれと表す。

コールのガンマ

デルタの計算結果を用いて以下のように計算することが出来る。

(1)式より、

プットのガンマ

(4)式より、

よって、

ベガ（ν）

コール/プットオプションのベガをそれぞれと表す。

コールのベガ

ここで、(1)式より

また、(3)式より

であるから

ここで、コールのガンマを求める際と同様に、第一項は0になる。従って

プットのベガ

(4)式より、

従って、

ロー（ρ）

コールのロー

ここで、

と(1)式より

コールのガンマを求める際と同様に、第二項は０であるので

プットのロー

(4)式より、

従って、

イプシロン（ε）

コールのイプシロン

ここで、

と(1)式より

コールのガンマを求める際と同様に、第二項は０であるので

プットのイプシロン

(4)式より

従って、

シータ（Θ）

コールのシータ

(1)式より

ここで、(3)式より

と(2)式より

プットのシータ

(4)式より

よって

従って

実装

pythonで各グリークスを実装する。

準備

共通する準備の部分

#グリークスを計算する。

#値の設定
K = 105;
r = 0.01;
q = 0.02;
sigma = 0.3;
tau = 0.1; #T - t

##d1,d2を計算する。
import math
from scipy.stats import norm
import matplotlib.pyplot as plt

def d1(S,K,r,q,sigma,tau):
    d1 = (math.log(S/K) + (r - q + sigma**2 / 2) * tau) / (sigma * math.sqrt(tau))
    
    return d1

def d2(S,K,r,q,sigma,tau):
    d1 = (math.log(S/K) - (r - q + sigma**2 / 2) * tau) / (sigma * math.sqrt(tau))
    
    return d1

#株価のリストを作成
S_list = list(range(80,120))
d1_list = []
d2_list = []

#d1を計算する。
for i in S_list:
    d1_list.append(d1(i,K,r,q,sigma,tau))

#d2を計算する。
for i in S_list:
    d2_list.append(d2(i,K,r,q,sigma,tau))

デルタ

#デルタを求める。
###コールのデルタ
def call_delta(q,T,d1):
    delta = math.exp(- q * T) * norm.cdf(d1, loc = 0, scale = 1)
    
    return delta
    
###プットのデルタ
def put_delta(q,T,d1):
    delta = - math.exp(- q * T) * norm.cdf(- d1, loc = 0, scale = 1)
    
    return delta


##コール、プットのデルタを計算し、格納する。
c_delta_list = []
p_delta_list = []

for i in range(len(S_list)):
    c_delta_list.append(call_delta(q,T,d1_list[i]))
    p_delta_list.append(put_delta(q,T,d1_list[i]))))

#デルタをプロットする。
import matplotlib.pyplot as plt
fig, ax = plt.subplots()

ax.plot(S_list, c_delta_list, label="CallDelta")
ax.plot(S_list, p_delta_list, label="PutDelta")
ax.set_title("Call and Put Delta")
ax.legend(loc = 0)

ax.set_xlabel("underlying asset price")
ax.set_ylabel("Delta")

plt.show();

ガンマ

#ガンマを求める
##ガンマを求める関数を定義
def gamma(S,sigma,tau,d):
    gamma = (1 / (S * sigma * math.sqrt(tau))) * (1 / math.sqrt(2 * math.pi)) * math.exp(- d ** 2 / 2) * math.exp(- q * tau)
    
    return gamma

##ガンマを実際に求め、値を格納する。
gamma_list = []
for i in range(len(S_list)):
    gamma_list.append(gamma(S_list[i],sigma,tau,d1_list[i]))

##ガンマをプロットする。
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
plt.plot(S_list, gamma_list, label="Gamma")
plt.title("Gamma")
plt.xlabel("underlying asset price")
plt.ylabel("Gamma")

plt.show()

ベガ

#ベガを計算する。
##ベガを計算する関数を定義する
def vega(S,sigma,tau,d1):
    vega = S * math.sqrt(tau / (2 * math.pi)) * math.exp( -d1 ** 2 / 2 - q * tau)
    
    return vega
##ベガを実際に求める。
vega_list = []

for i in range(len(S_list)):
    vega_list.append(vega(S_list[i],sigma,tau,d1_list[i]))

##ベガをプロットする。
fig, ax = plt.subplots()
plt.plot(S_list, vega_list, label="Vega")
plt.title("Vega")
plt.xlabel("underlying asset price")
plt.ylabel("Vega")

plt.show()

ロー

#ローを計算する。
##ローを計算する関数を定義する

###コールのロー
def call_rho(K,r,tau,d2):
    rho = K * tau * math.exp(-r * tau) * norm.cdf(d2, loc = 0, scale = 1)
    
    return rho

###プットのロー
def put_rho(K,r,tau,d2):
    rho = -K * tau * math.exp(-r * tau) * norm.cdf(-d2, loc = 0, scale = 1)
    
    return rho


##ローを実際に求める。
c_rho_list = []
p_rho_list = []

for i in range(len(S_list)):
    c_rho_list.append(call_rho(K,r,tau,d2_list[i]))
    p_rho_list.append(put_rho(K,r,tau,d2_list[i]))

#ローをプロットする。
fig, ax = plt.subplots()

ax.plot(S_list, c_rho_list, label="CallRho")
ax.plot(S_list, p_rho_list, label="PutRho")
ax.set_title("Call and Put Rho")
ax.legend(loc = 0)

ax.set_xlabel("underlying asset price")
ax.set_ylabel("Rho")

plt.show();

イプシロン

#ローエフを計算する。
##ローエフを計算する関数を定義する

###コールのローエフ
def call_rhof(S,q,tau,d1):
    rhof = - S * tau * math.exp(- q * tau) * norm.cdf(d1, loc = 0, scale = 1)
    
    return rhof

###プットのローエフ
def put_rhof(S,q,tau,d1):
    rhof = S * tau * math.exp(- q * tau) * norm.cdf(- d1, loc = 0, scale = 1)
    
    return rhof


##ローエフを実際に求める。
c_rhof_list = []
p_rhof_list = []

for i in range(len(S_list)):
    c_rhof_list.append(call_rhof(S_list[i],q,tau,d1_list[i]))
    p_rhof_list.append(put_rhof(S_list[i],q,tau,d1_list[i]))

#ローエフをプロットする。
fig, ax = plt.subplots()

ax.plot(S_list, c_rhof_list, label="CallRhof")
ax.plot(S_list, p_rhof_list, label="PutRhof")
ax.set_title("Call and Put Rhof")
ax.legend(loc = 0)

ax.set_xlabel("underlying asset price")
ax.set_ylabel("Rhof")

plt.show();

シータ

#シータを計算する。

##コールのシータ
def call_theta(S,K,r,q,sigma,tau,d1,d2):
    theta = q * S * math.exp(- q * tau) * norm.cdf(d1, loc = 0, scale = 1) - r * K * math.exp(- r * tau) * norm.cdf(d2, loc = 0, scale = 1) - S * sigma * math.exp(-q * tau) / (2 * math.sqrt(tau)) * (1 / math.sqrt(2 * math.pi)) * math.exp( - d1 ** 2 / 2)
    
    return theta

###プットのシータ
def put_theta(S,K,r,q,sigma,tau,d1,d2):
    theta = - q * S * math.exp(- q * tau) * norm.cdf(- d1, loc = 0, scale = 1) + r * K * math.exp(- r * tau) * norm.cdf(- d2, loc = 0, scale = 1) - S * sigma * math.exp(-q * tau) / (2 * math.sqrt(tau)) * (1 / math.sqrt(2 * math.pi)) * math.exp( - d1 ** 2 / 2)
    
    return theta


##シータを実際に求める。
c_theta_list = []
p_theta_list = []

for i in range(len(S_list)):
    c_theta_list.append(call_theta(S_list[i],K,r,q,sigma,tau,d1_list[i],d2_list[i]))
    p_theta_list.append(put_theta(S_list[i],K,r,q,sigma,tau,d1_list[i],d2_list[i]))    


#シータをプロットする。
fig, ax = plt.subplots()

ax.plot(S_list, c_theta_list, label="CallTheta")
ax.plot(S_list, p_theta_list, label="PutTheta")
ax.set_title("Call and Put Theta")
ax.legend(loc = 0)

ax.set_xlabel("underlying asset price")
ax.set_ylabel("Theta")

plt.show();

何か間違い等あればコメントで教えてください。

おわり

---

1. 一般的には配当率qに対する感応度に対応するグリークスはや、で知られており、ここではとする。ちなみに大学院では"Rhof"（どのギリシャ文字に対応するのかは不明）で習った。↩

コイン投げから一様分布（と正規分布）を生成する方法

表題の通り
ファイナンスのための確率解析Ⅱ（S.E.シュリーブ）の第一章に書いてあった一様分布の生成方法を検証してみる。

具体的な方法

表が出る確率が1/2であるコインを考える。（要はp = 1/2のベルヌーイ試行）
このコインをn回投げるとする。n = 1,2,...に対して
\[
Y_n(\omega) = \begin{cases}
1 & \omega_n = H \\
0 & \omega_n = T
\end{cases}
\]
そして、
\[
X = \sum^{\infty}_{n=1} \frac{Y_n}{2^n}
\]
と定義する。

この時、Xは[0,1]上に一様に分布する。

検証

MATLABを用いて一様分布に従っているか調べる。残念ながら私は高々有限回しかコインを投げられないので*1、とりあえずn = 5,10,100と値を変えて調べることにする。また、サンプルサイズは1000とする。

結果

階級の幅を0.1としたヒストグラムは以下の通りである。

また、縦軸に経験分布、横軸を理想的な一様分布としたqqプロット*2は以下の通りである。

qqプロット

プロットされたデータを見ると、n=5でも十分に一様分布しているように見える。（一方でnを増やすと、より一様分布するとは限らない？）

また、定性的な結果を得るために上のヒストグラムを作成した際に離散化した乱数値に対してカイ二乗適合度検定*3を行う。自由度は９である。
有意水準5％において、すべてのnに対して帰無仮説は棄却されなかった。従って、離散化した乱数値は離散一様分布に従っていないとは言えない*4ということが分かった。*5

更なる検証

5回コインを投げればそれなりの一様乱数を生成できることが分かったので、これを用いて標準正規乱数を作成してみる。ここでは、ボックス＝ミュラー法*6を用いることにする。
尚、ボックス＝ミュラー法では一様乱数が二つ必要なので、乱数は二倍必要とされる。（n=5であれば10回、n=10であれば20回コインを投げる必要がある。）

結果は以下の通りである。

ヒストグラムを見る限りでは正規分布に従っているようだが、qqプロットによると、裾の方が標準正規分布から乖離していることがわかる。

また、生成された乱数に対してコルモゴロフースルミノフ検定*7を行ったところ、帰無仮説は棄却された。従って、標準正規乱数であるとは言えないことが分かった。

nの値を増やしてもう一度検証をしてみる。とりあえずn=10とする。

今度はqqプロットの結果も含め正規分布に従っているように見える。実際、コルモゴロフ-スルミノフ検定の結果は棄却されなかった。

まとめ

コイン投げの結果から一様乱数、更には標準正規乱数を生成してみようと試みた。結果としては、単に一様分布を生成したいだけであればn=5でもある程度十分な精度の乱数が得られる事が分かった。
ただ、そこから更に標準正規乱数を生成しようとするとn=5では不十分であり、n=10、つまり20回コインを投げる必要がある事が分かった。(20回もコイン投げるんだったら二項分布から正規分布への近似法使った方がよくね？？）

もし無人島に裸一貫で飛ばされたとしても、（歪んでいない）コインが一枚あれば一様分布、正規分布、さらには逆関数法を使えば他にも様々な分布に従う乱数が作れると思うと夢が膨らみますね！！（強引な結論）

使用したコード（matlab）

%コイン投げから一様分布を生成し、プロットや検定を行うコード
n_list = [5, 10, 100] %ループ回数のリスト

for i = 0:2
    
    n = n_list(1, i + 1) %ループ回数を設定
    
    sample_size = 1000;
    
    %乱数の配列を作る（sample_size × 1)
    ransu = double.empty(sample_size,0); %初期配列
    for j = 1:sample_size
        ransu(j,1) = itiyo(n);
    end
    
　%ヒストグラムを作成する
    figure
    histogram(ransu)
    xlabel("乱数の値")
    ylabel("度数")
    title(['乱数の分布 n = ', num2str(n)], 'FontSize',12)
    
    %乱数が一様分布に従っているか、カイ二乗適合度検定を行う。(自由度は9とする。）
    %%カイ二乗統計量を求める
    toukeiryo = kai2_kentei(sample_size, ransu)
    
    %%上側確率を求める。
    upper_prob = 1 - chi2cdf(toukeiryo,9)
    
    %%結果をansに入れる。
    if upper_prob >= 0.05
        ans = "一様分布していないとはいえない" %棄却できない
    else
        ans = "一様分布していない" %棄却できる
    end
    
    %qqプロットを行う
    pd = makedist('Uniform');
    
    figure;
    qqplot(ransu,pd);
    title(['標本データvs一様分布（理想）のqqプロット n = ', num2str(n)])
end

%以下関数
function itiyo = itiyo(n) %ベルヌーイ分布から一様分布を作成する関数、nはループ回数
    sum = 0;
    
    for i = 1:n
        ber = binornd(1, 0.5);
        sum = sum + ber / (2^i);    
    end

    itiyo = sum;
end


function kai2_kentei = kai2_kentei(sample_size, ransu) %一様分布に従っているか、カイ二乗適合度検定のための統計量を算出する関数
    %初期値
    kai = 0;
    
    for i = 0:9
        expect = sample_size * 0.1; %期待度数
        obs = length(ransu( 0.1 * i < ransu & ransu < 0.1 * i + 0.1)); %観測度数
        kai = kai + (expect - obs)^2 / expect; %カイ二乗統計量を求める。
    end
    
    kai2_kentei = kai;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%box-muller法を使用し、標準正規乱数を生成するコード

n = 10

sample_size = 1000;

%一様乱数を生成する(samplesize ×2の配列）
for i = 1:2
    for j = 1:sample_size
        ransu(j,i) = itiyo(n);
    end
end

%box mullar法を用いて、一様乱数から正規乱数を生成する。
for s = 1 : sample_size
    X1(s) = sqrt( -2 * log(ransu(s,1))) * cos(2 * pi * ransu(s,2));
end

histogram(X1)
xlabel("乱数の値")
ylabel("度数")
title(['乱数の分布 n = ', num2str(n)], 'FontSize',12)

%qqプロットを行う
pd = makedist('Normal');

figure;
qqplot(X1,pd);

%コルモゴロフ-スミルノフ検定
h = kstest(X1)

if h == 0
    ans = "正規分布に従っていないとは言えない"
else
    ans = "正規分布に従っていない"
end

%以下関数
function itiyo = itiyo(n) %ベルヌーイ分布から一様分布を作成する関数、nはループ回数
    sum = 0;
    
    for i = 1:n
        ber = binornd(1, 0.5);
        sum = sum + ber / (2^i);    
    end

    itiyo = sum;
end

*1:可算無限回コインを投げられる人は各自で実験してみて下さい。

*2:もし乱数が一様分布に従っているのであれば、プロットされた点はほぼ45度線上に分布する。詳しくは以下のURLを参照 ja.wikipedia.org

*3:帰無仮説は「度数が（離散）一様分布に従っている」

*4:ぺこぱ

*5:あくまで離散化したデータに対してそう言えるだけで、元の連続的な乱数値に対しては何も言えないのに注意（階級の幅の取り方によって結果は変わる）

*6: (0,1)上の一様分布に従う確率変数X,Yを用いて、 \begin{align} Z_1 =\sqrt{-2 logX} cos 2 \pi Y \\ Z_2 =\sqrt{-2 logX} sin 2 \pi Y \\ \end{align} とすれば、Z_1とZ_2は独立に標準正規分布に従うというもの。 詳しくは ja.wikipedia.org

*7:ベクトルデータが標準正規分布に従っているか調べる検定。帰無仮説は「データが標準正規分布から発生している。」 jp.mathworks.com

対数ロジスティック分布

はじめに

統計検定一級（理工学）の為、「一般化線形モデル入門　原著第２版（Annette J. Dobson著)」を使って生存時間解析の勉強をしていたら、「対数ロジスティック分布」という聞いた事があるようでよく知らない確率分布が出現。
ネットで調べてみるも日本語のサイトが全くと言っていいほど見つからず、結局英語版WikipediaやらMATLAB公式サイトの記事等を使って手計算するしかなかった。
以下、備忘録として途中式を中心に記録を残しておきます。

対数ロジスティック分布（log-logistic distribution）とは

対数を取ると「ロジスティック分布」に従うような分布を「対数ロジスティック分布」という。確率変数が非負の値を取り、かつ密度関数の裾が対数正規分布よりも厚いことから、生存時間解析の分野で使われることが多い（らしい）

ロジスティック分布（logistic distribution）

正規分布と同じく釣鐘型の密度関数を持つ連続型確率分布。正規分布よりも裾が厚い（平均から離れても下がりにくい）。

パラメータ(μ,s)はそれぞれ位置 (location) パラメータと尺度(scale)パラメータと呼ばれ、期待値、分散はそれぞれ
\begin{align}
E[X] &= \mu \\
V[Y] &= \frac{\pi ^2 s^2}{3}
\end{align}
である。また、累積分布関数は
\begin{equation}
F_X(x) = \frac{1}{1 + e ^ {-\frac{(x-\mu)}{s}}}
\end{equation}
である。

以下のグラフは標準正規分布と標準ロジスティック分布（(μ,s) = (0,1)となるようなロジスティック分布）の確率密度関数を比較したものである。

累積分布関数(cumulative distribution function)

上で挙げたロジスティック分布の累積分布関数から計算することが出来る。
Xがロジスティック分布に従い、Yが対数ロジスティック分布に従うとすると
$$
log(Y) \overset{\mathrm{dist}}{=} X \Leftrightarrow Y \overset{\mathrm{dist}}{=} exp(X)
$$
であることから、Yの累積分布関数は
\begin{align}
F_Y(y) &= P(Y \leq y) \\
&= P(e^X \leq y) \\
&= P(X \leq log(y)) \\
&= F_X(log y) \\
&= \frac{1}{1 + exp\{ -\frac{(logy - \mu)}{s}\}} \\
&= \frac{1}{1 + exp\{(logy^{-\frac{1}{s}} + \frac{\mu}{s}\}} \\
&= \frac{1}{1 + y^{-\frac{1}{s}} e^{\frac{\mu}{s}}} \\
&= \frac{1}{1 + \left( \frac{e^{\mu}}{y} \right)^{\frac{1}{s}}} \\
&= \frac{1}{1 + \left( \frac{y}{e^{\mu}} \right)^{-\frac{1}{s}}} \tag{1}　
\end{align}

となる。このままでも良いが、WikipediaやTadikamalla, Pandu R. (1980)では
$$\mu = log(\alpha), s = \frac{1}{\beta}$$
としてパラメータを変換した形で書いている。つまり累積分布関数は
\begin{align}
F_Y(y) &= \frac{1}{1 + \left( \frac{y}{\alpha} \right)^{-\beta} } \tag{2} \\
&= \frac{y^{\beta}}{y^{\beta} + \alpha^{\beta}} \tag{3}
\end{align}
と表すことが出来る。以下、この形式を定式として進める。

確率密度関数(probability density function)

上で得られた累積分布関数(2)式をyで微分すればよい。つまり、

\begin{align}
f_Y(y) &= - \frac{ \left( 1 + \left( \frac{y}{\alpha} \right)^{-\beta} \right)^{'} } {\left( 1 + \left( \frac{y}{\alpha} \right)^{-\beta} \right)^2 } \\
&= \frac{ \left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{-\beta - 1}} {\left( 1 + \left( \frac{y}{\alpha} \right)^{-\beta} \right)^2 } \\
&= \frac{ \left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{-\beta - 1} \left( \frac{y}{\alpha} \right) ^ {2 \beta}} {\left( 1 + \left( \frac{y}{\alpha} \right)^{-\beta} \right)^2 \left( \frac{y}{\alpha} \right) ^ {2 \beta}} \\
&= \frac{\left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{\beta - 1}} {\left( 1 + \left( \frac{y}{\alpha} \right)^{\beta} \right)^2 } \tag{4}　
\end{align}

生存関数（survivor function）

生存時間解析における「生存関数（survivor function）」とは、Yを生存時間とした時、Yがある時間yを超える確率を意味している。つまり、

\begin{align}
S_Y(y) &= P(Y \geq y)\\
&= 1 - F_Y(y)
\end{align}
と定義される。

従って、対数ロジスティック分布の生存関数は(3)式より
\begin{align}
S_Y(y) &= 1 - \frac{y^{\beta}}{y^{\beta} + \alpha^{\beta}} \\
&= \frac{\alpha^{\beta}}{y^{\beta} + \alpha^{\beta}} \\
&= \frac{1}{1 + \left( \frac{y}{\alpha} \right) ^{\beta}} \tag{5}
\end{align}
となる。

ハザード関数（hazard function）

生存時間解析における「ハザード関数（hazard function）」とは、Yを生存時間とした時、時間yまでは生存しているという条件のもとでyとy + δyの間に死亡する確率の極限値（δy→０）を意味している。
すなわち、
\begin{align}
h(y) &= \frac{f(y)}{S(y)} \tag{6} \\
\end{align}
と定義される。また、$$f(y) = -\frac{dS(y)}{dy}$$であることから、
\begin{align}
h(y) &= -\frac{S(y)^{'}}{S(y)} \\
&= -\frac{d}{dy} (log[S(y)]) \tag{7}
\end{align}
という表し方も出来る。

対数ロジスティック分布のハザード関数は、(6)(5)(4)式より
\begin{align}
h_Y(y) &= \frac{f_Y(y)}{S_Y(y)} \\
&= \dfrac{\frac{\left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{\beta - 1}} {\left( 1 + \left( \frac{y}{\alpha} \right)^{\beta} \right)^2 }}{\frac{1}{1 + \left( \frac{y}{\alpha} \right) ^{\beta}}} \\
&= \frac{\left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{\beta - 1}} {1 + \left( \frac{y}{\alpha} \right)^{\beta}}
\end{align}
となる。

k次モーメント

k次のモーメントは、上で求めた確率密度関数を用いて以下のように計算できる。

\begin{align}
E[Y^k] &= \int^\infty_0 y^k f(y) dy \\
&= \int^\infty_0 y^k \frac{\left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{\beta - 1}} {\left( 1 + \left( \frac{y}{\alpha} \right)^{\beta} \right)^2 } dy \\
\end{align}

$$ここで、t = \left( \frac{y}{\alpha} \right)^{\beta} とすると、$$
$$dt = \left( \frac{\beta}{\alpha} \right) \left( \frac{y}{\alpha} \right)^{\beta -1} dy$$

また自然数kを用いて
\begin{align}
t^{\frac{1}{\beta}} &= \frac{y}{\alpha}より、\\
y^k &= (\alpha t^{\frac{1}{\beta}})^k
\end{align}
であることから、

\begin{align}
E[Y^k] &= \int^\infty_0 (\alpha t^{\frac{1}{\beta}})^k \frac{1}{(t + 1)^2} dt \\
&= \alpha^k \int^\infty_0 \frac{t^{\frac{k}{\beta}}}{(t + 1)^2} dt \\
&= \alpha^k \int^\infty_0 \frac{t^{(1 + \frac{k}{\beta}) - 1}}{(t + 1)^{(1 + \frac{k}{\beta})+(1 - \frac{k}{\beta})}} dt \\
&= \alpha^k B(1 - \frac{k}{\beta}, 1 + \frac{k}{\beta}) \tag{8}
\end{align}
が得られる。ただし、$$B(・,・)$$はベータ関数。
（Wikipediaだとここから更に
$$ E[Y^k] = \alpha^k \frac{\frac{k \pi}{\beta}}{sin(\frac{k \pi}{\beta})}$$
という風に変形していますが仕組みがよく分かっていません...　分かる方は是非コメントで教えて下さい）

これにより、平均、分散はそれぞれ
\begin{align}
E[Y] &= \alpha B(1 - \frac{1}{\beta}, 1 + \frac{1}{\beta})\\
V(Y) &= E[Y^2] - E[Y] \\
&= \alpha^2 B(1 - \frac{2}{\beta}, 1 + \frac{2}{\beta}) - \alpha B(1 - \frac{1}{\beta}, 1 + \frac{1}{\beta}) \\
\end{align}
で求めることが出来る。

余談

各グラフはMATLABで描いたのだが、MATLABにおける対数ロジスティック分布は(1)式の形で定義されているので、(2)式の形で手計算した場合パラメータ変換の注意が必要である。

参考文献

・Annette J. Dobson 著、田中豊、森川敏彦、山中竹春、冨田誠 訳「一般化線形モデル入門 原著第２版」
https://www.amazon.co.jp/%E4%B8%80%E8%88%AC%E5%8C%96%E7%B7%9A%E5%BD%A2%E3%83%A2%E3%83%87%E3%83%AB%E5%85%A5%E9%96%80-%E5%8E%9F%E8%91%97%E7%AC%AC2%E7%89%88-Annette-J-Dobson/dp/4320018672
・Tadikamalla, Pandu R. (1980), A Look at the Burr and Related Distributions, International Statistical Review, 48 (3)
・Log-logistic distribution - Wikipedia
・対数ロジスティック分布 - MATLAB & Simulink - MathWorks 日本