初投稿です。
至らない点等あるかもしれませんがよろしくお願いします。

概要

Juliaを用いて、株価の変動がブラック・ショールズモデルに従うとしたときのシミュレーションを行います。
本記事においては、~~筆者に厳密な議論をするだけ数学力が無いので~~種々の厳密さを抜きにして"お気持ち"だけで進んで行きます。

対象読者

正規分布が(ふんわり)わかる人
Juliaの基本的な文法(forとかifとか)が一通りわかるけど、試しになにしよう?ってなってる人

環境

Julia 1.5.3
使用パッケージ
- Distributions 0.24.4
- Plots 1.20.1

株価の変動モデル

時刻$t$における株価を$S_t$とします。$S_t$の変動が以下のような幾何ブラウン運動に従うとしたのがブラック・ショールズモデルです。

\begin{align}
dS_t = a S_t dt + \sigma S_t dW_t \tag{1}
\end{align}

ここで、

$a$は期待収益率(定数) ⇒平均的にどれだけ収益率があるか
$\sigma$はボラティリティ(定数) ⇒株価の変動がどれくらい激しいか
$W_t$はブラウン運動(後述)

です。
(1)式のお気持ちを説明すると、株価$S_t$の変動$dS_t$を

全体的な傾向を表す部分(第1項)
確率的な変動を表す部分(第2項)

の2つの部分の和として表現している感じですね。

さて、$W_t$とおいたブラウン運動ですが、それは以下のような性質を持つ確率過程です。

出発位置
$W_0 = 0$
連続性
$W_t$はほとんど確実に連続。
正規増分
任意の$0 \leq s < t$に対して、$W_t - W_s$が平均0、分散$t-s$の正規分布に従う。
独立増分
任意の$0=t_1 \leq t_2 < t_3 < \dots < t_n $に対して
$W_{t_2} - W_{t_1} , W_{t_3} - W_{t_2} , \dots , W_{t_n} - W_{t_{n-1}}$は独立。

シミュレーションをする上で押さえておきたいのは、後半の2つ(正規増分と独立増分)です。
後半の二つを~~都合良く解釈する~~考えると、$dW_t$については以下のように考えてシミュレーションに落とし込めば良さそうです(本当か?)
$$dW_tはN(0,dt)に従う$$

対数価格と伊藤の補題

さて、今までの考察から既に株価$S_t$の変動をシミュレーションすることも可能なのですが、金融工学等の世界では(1)式のまま株価を扱うことはあまりありません。
これから紹介するように、対数価格$X_t$を興味の対象にすることが多い印象です。
~~???「それあなたの感想ですよね?」~~

なぜ価格の対数を取るのか

色々理由はありますが、個人的に納得した理由は

対数価格の差が収益率の近似になっているから
対数価格の変動はすっきりしてるから

の2点ですね。
まず、１点目について説明しましょう。
収益率$R_{t_i}$は
$$
R_{t_i} = \frac{S_{t_i}-S_{t_{i-1}}}{S_{t_{i-1}}}
$$
と表されます。ここで$\log{(1+x)}$を$x=0$周りでテイラー展開して、一次の項までで近似すると
$$
\log{(1+x)} \simeq x
$$
となります。この$x$に$x = R_{t_i}$を代入すると、

R_{t_i} \simeq \log{ \Big( 1+\frac{S_{t_i}-S_{t_{i-1}}}{S_{t_{i-1}}} \Big)} \\
= \log{S_{t_i}} - \log{S_{t_{i-1}}} \\
\therefore R_{t_i} \simeq X_{t_i} - X_{t_{i-1}}

となり、収益率の近似になっていることが確認できます。

次に2点目なのですが、これを説明するためには伊藤の補題について説明する必要があります。

伊藤の補題

伊藤の補題(或いは伊藤の公式)とは、以下のようなものです。

$S_t$は下記の確率過程(伊藤過程)に従うとする。
$$
dS_t = a_t dt + b_t dW_t \tag{2}
$$
また、$t,x$の関数$f(t,x)$が$t$について1回微分、$x$について2回微分可能とする。
$f(t,x)$の偏導関数を

f_t = \frac{\partial f}{\partial t} \\
f_x = \frac{\partial f}{\partial x} \\
f_{xx} = \frac{\partial^2 f}{\partial x^2}

とする。
この時、

\begin{align}
df(t,S_t) =& \Big( f_t(t,S_t) + f_x(t,S_t)a_t + \frac{1}{2} f_{xx} (t,S_t) b^2_t \Big)dt \\
&+ f_x(t,S_t)b^2_t dW_t \tag{3}
\end{align}

となる。

さて、これを使って(1)式をもとに対数価格$X_t$の変動を考えてみましょう。
$f(t,x) = \log{x}$とすると、

\begin{align}
f_t(t,x)= \frac{\partial f}{\partial t} &= 0 \\
f_x(t,x) = \frac{\partial f}{\partial x} &= \frac{1}{x} \\
f_{xx}(t,x) = \frac{\partial^2 f}{\partial x^2} &= -\frac{1}{x^2}
\end{align}

であり、(1)式と(2)式の対応関係は

(2)式の$a_t$:$a S_t$
(2)式の$b_t$:$\sigma S_t$

なので、(3)式より

df(t,S_t) = dX_t = \Big( a - \frac{\sigma^2}{2} \Big)dt + \sigma dW_t

となります。ここで$\mu \equiv \Big( a - \frac{\sigma^2}{2} \Big)$とおけば、

$$
dX_t = \mu dt +\sigma dW_t \tag{4}
$$

となって対数価格の変動がすっきりします。

数値シミュレーション

お待たせしました。
今までの議論妄想の垂れ流しを用いて、対数価格$X_t$をシミュレーションして図示してみましょう。

時間の分割設定

時間の分割を下記のように等間隔に定めます。
(Juliaは配列のIndexが1始まりなので、時間のIndexも1から始めておきましょう。)

\begin{align}
&0 = t_1 < t_2 < \dots < t_n = T \\
&t_i = \frac{T}{n-1} (i-1) \quad (i \in \{1, 2, \dots, n \})\\
&\Big( \Delta t = \frac{T}{n-1} \Big) \\
&(但し n \geq 2 かつ T > 0)
\end{align}

ただ等間隔に区切っただけですね。これで時間の設定ができました。
いよいよコードに入っていきます。

コード

まず、コードの全体像をお見せします。

logPricePlot.jl

using Plots
using Random
using Distributions

function logPricePlot(n,seed = 10)
    #---時間の生成---
    T = 1
    Δt = T/(n-1)
    t = collect(0:Δt:T) #一定間隔Δtで、Tまで時間を生成

    #---価格の生成---
    #正規分布を生成
    d = Normal(0,Δt^0.5) #dWは平均0,標準偏差Δt^0.5に従う

    #正規分布に従う乱数を生成
    #Random.seed!(seed) #乱数のシードを固定したい人向け
    ΔW = rand(d,n-1)

    #対数価格変動をシミュレート
    μ = 0.04
    σ = 0.2
    X = zeros(n)
    X[1] = 10 #初期対数価格X_1を設定
    for i in 1:n-1
       X[i+1] = X[i] + μ*Δt + σ*ΔW[i] 
    end

    #---価格のプロット---
    plot(t,X,
        xlabel = "t",#x軸のラベル設定
        ylabel = "logPrice",#y軸のラベル設定
        label = "logPrice",#凡例設定
    )

    #グラフを保存したければ
    #savefig("画像名.png")
end

logPricePlot(1000) #データ数1000個でプロット

これを実行すると以下の感じに対数価格が生成できていると思います。

順に釈明説明していきたいと思います。

時間の生成

$T=1$と定めて時間を等間隔に生成しています。
このパートで分かりにくいのは

t = collect(0:Δt:T)

の部分だと思います。これは0から始めて間隔$\Delta t$で$T$までのデータを作って、collect()で配列に直しているという感じですね。
配列にしなくても別に大丈夫で、単に

t = 0:Δt:T

でも動きます。

価格の生成

ここがメインパートです。
どういうことを行っているかというと(4)式
$$
dX_t = \mu dt +\sigma dW_t \tag{4}
$$
を以下のように読み替えて、逐次的に対数価格を作っています。

\begin{align}
dX_t &= \mu dt +\sigma dW_t \tag{4} \\
\Delta X_{t_i} &= \mu \Delta t + \sigma \Delta W_{t_i} \\
X_{t_{i+1}} - X_{t_i} &= \mu \Delta t + \sigma \Delta W_{t_i} \\
X_{t_{i+1}} &= X_{t_i} + \mu \Delta t + \sigma \Delta W_{t_i} \\
(\Delta W_{t_i} &\sim N(0,\sqrt{\Delta t}))
\end{align}

$\Delta W_{t_i}$を作る部分、即ち正規分布に従う乱数を生成する部分では、Juliaの標準パッケージであるRandomと、外部パッケージであるDistributionsの力を借りています。

正規分布に従う乱数が欲しい場合には、まず正規分布を作ります。

d = Normal(0,Δt^0.5) #平均0,標準偏差Δt^0.5の正規分布を作成

次に作成した分布をrand()に渡すことで、その分布に従う乱数を指定した個数作ります。

ΔW = rand(d,n-1) #rand(分布,(次元))で分布に従う乱数を指定した次元で作ってくれる

こうして$\Delta W_{t_i}$を作成したら、あとは逐次的に対数価格を作ってしまえば良いです。

価格のプロット

これはPlotsというパッケージの力を借りて、時間と対数価格をプロットしています。
特に解説が必要な部分とかはないと思います。

終わりに

今回はJuliaで対数価格のシミュレーションを行ってみました。
次回は~~書く気力があれば~~今回のシミュレーションをもとにRealized VolatilityをJuliaで計算してみたとかを書いてみようかと思ってます。

Juliaで株価のシミュレーション

概要

目次