금융공학에 대해 친절하게 잘 설명해 놓은 블로그 (Finance Diary)를 보고 Julia 버전으로 만들어 보았습니다.

SymPy를 좀 더 알겸 해서 SymPy로 코딩을 해봤습니다.

stock_math_model

주가의 수학적 모델¶

In [1]:

import SymPy as sp
import PyCall
using LaTeXStrings
using Latexify
import Markdown as mk

using StatsBase
using Distributions
using StatsPlots

using Random
using Dates
# sympy.stats를 sp에 import
stats = PyCall.pyimport_conda("sympy.stats","sympy")
sp.import_from(stats)

주요 참고¶

SymPy 참고¶

sympy.stats 참고¶

https://github.com/JuliaPy/SymPy.jl/issues/262

버전정보¶

In [2]:

versioninfo()

Julia Version 1.8.3
Commit 0434deb161e (2022-11-14 20:14 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 28 × Intel(R) Xeon(R) CPU E5-2690 v4 @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, broadwell)
  Threads: 16 on 28 virtual cores
Environment:
  LD_LIBRARY_PATH = /usr/local/cuda/extras/CUPTI/lib64:/home/shpark/julia-1.8.3/lib:/home/shpark/julia-1.8.3/lib/julia:/home/shpark/gurobi752/linux64/lib:/lib:/home/shpark/artelys/knitro-12.1.1-Linux-64/lib:/usr/local/cuda-11.1/lib64::/usr/local/cuda/lib64/:/usr/local/cuda/extras/CUPTI/lib64
  JULIA_HOME = /home/shpark/julia-1.8.3
  JULIA_NUM_THREADS = 16

LaTex 참고¶

LaTex : Equation numbering, reference

\label : name describing he equation
\tag : the label appearing next to the equationcan be a number or letters
\eqref : reference to the labeled equation
예) \label{eq:stock_return} \tag{1}

LeTex Symbol\ http://www-ph.postech.ac.kr/~bimin/Latex.htm \ https://bekaykang.github.io/posts/markdown-latex/

Rate of Return Modeling¶

참고 : https://sine-qua-none.tistory.com/27
$\mu$ : 1년 기대 수익율(평균 수익률)
$\sigma$ : 1년 수익률의 표준편차 (변동성)
$\Delta t$ : 극히 짧은 시간 간격 예) $\Delta t$가 $\frac{1}{365}$ 이라면 $\Delta t$는 1일 이라는 뜻이고 1일 변동성 = $\sigma \sqrt{\Delta t}$
$S_t$ : $t$ 시간에 주식 가격
$\Delta t$ 시간동안 주식 수익률 ($R_{\Delta t}$) $$ \frac{S_{t+\Delta t}-S_t}{S_t}=\frac{\Delta S_t}{S_t}, \;\; \Delta S_t = S_{t+\Delta t} – S_t\\ $$
$\Delta t$ 시간동안 주식의 표준편차 $$\sigma \sqrt{\Delta t}$$
$\Delta t$ 시간동안 주식의 기대수익률 $$ \mu \sqrt{\Delta t} $$
종합 하면 $R_{\Delta t}$는 평균이 $\mu \Delta t$이고 분산이 $\sigma \sqrt{\Delta t}$인 정규분포를 따른다. $$ R_{\Delta t} = \frac{S_{t+\Delta t} – S_t}{S_t} \sim N \left(\mu \Delta t, \sigma \sqrt{\Delta t} \right) \label{eq:rate_of_return_1} \tag{1} $$
$Z_t$를 정규분포 $N \left(0,\Delta t \right)$을 따르는 확률변수라고 하면 $$ R_{\Delta t} = \frac{\Delta S_t}{S_t} = \mu \Delta t + \sigma Z_t, \;\;Z_t \sim N \left(0,\Delta t \right) \label{eq:rate_of_return_2} \tag{2} $$

$\Delta t$와 $\Delta S_t$를 무한히 작은 값(미분소)으로 바꾸면 $$ \frac{dS_t}{S_t} = \mu dt + \sigma Z_t, \;\; Z_t \sim N(0,dt) \label{eq:rate_of_return_3} \tag{3} $$

위너 프로세스(Wiener process)¶

확률과정(Stochastic process ) : 시간에 따라서도 계속정의 되는 확률변수 ### 위너 프로세스의 4가지 성질

임의의 시점 $t > 0$와 time interval $\Delta t > 0$에 대하여, $W_{t+\Delta t} – W_t$는 $N\left(0,\Delta t \right)$을 따른다.
(현재)시점 $t=0$에서 이 값은 0 이다. 즉, $W_0 = 0$ 이다.
임의의 시점 $t_1 \; < t_2 \; < \; t_3$ 에 대하여 $W_{t3} – W_{t2}$ 와 $W_{t1}$ 은 독립 이다.\ (즉, 시점 $t_1$ 전까지 발생한 사건들은 그 이후 시점에 벌어진 사건들과 관계가 없다.)
$t$에 대해서 연속함수 이다.

위너 프로세스를 사용하여 수식 $\eqref{eq:rate_of_return_3}$을 다시 모델링 하면\ 위너 프로세스의 1번 성질을 통해 $$ Z_t \sim N\left(0,dt\right)\\ W_{t+dt} – W_t \sim N\left(0,dt\right)\\ dW_t = W_{t+dt} – W_t = Z_t \sim N\left(0,dt\right)\\ \\ \frac{dS_t}{S_t} = \mu dt + \sigma dW_t $$

Expectation¶

In [3]:

sp.@vars x X μ Z dt dW t W S dS S₀ σ T

Out[3]:

(x, X, μ, Z, dt, dW, t, W, S, dS, S₀, σ, T)

In [4]:

pdf(x::sp.Sym,s::sp.Sym) = sp.density(x).pdf(s)

Out[4]:

pdf (generic function with 1 method)

In [5]:

pdf(sp.Normal(x,μ,σ),x)

Out[5]:

$\frac{\sqrt{2} e^{- \frac{\left(x – μ\right)^{2}}{2 σ^{2}}}}{2 \sqrt{\pi} σ}$

In [6]:

X = sp.Normal(x,μ,σ)
sp.E(X^4)(μ=>0)

Out[6]:

$3 σ^{4}$

$dW_t \sim N(0,\sqrt{dt})$
$x = dW_t$
$\sigma = \sigma_1 \sqrt{dt}$
$\mu=0$
$x$ 가 정규분포를 따른다고 할 때 $x^2$의 분산 \ $Var(x^2)$ = $E\left((x^2)^2\right) – \left(E(x^2)\right)^2$ = $E\left(x^4\right) – \left(E(x^2)\right)^2 = 2dt^2$

In [7]:

dW = sp.Normal(x,0,√dt)
"dW² 의평균" |> display
sp.E(dW^2) |> display
"dW² 의 분산 " |> display
v = sp.variance(dW^2) 
v |> display

"dW²의 분산 var 가 dt 가 0으로 가는 경우 0 이 되기 때문에 dW² 자체가 dt 가 된다" |> display
sp.limit(v,dt, 0) |> display
"따라서 dW²= dt" |> display

"dW² 의평균"

$dt$

"dW² 의 분산 "

$2 dt^{2}$

"dW²의 분산 var 가 dt 가 0으로 가는 경우 0 이 되기 때문에 dW² 자체가 dt 가 된다"

$0$

"따라서 dW²= dt"

Taylor Series¶

In [8]:

let
    sp.@vars S x h dS x0 ΔS Δt ΔW     
    sp.@symfuns f
    F = f(S+ΔS).series(ΔS,n=3)
    F |> display
    F = F.removeO()
    F |> display
    F = F.doit()
    F |> display
    F = F(ΔS=>(μ*Δt + σ*ΔW)*S).simplify()
    F |> display
#     sp.simplify(F(S=>log(S))) |> display
    F = F-f(S)
    F = F.doit()
    F |> display
    L"f(S+ΔS)-f(S)=%$(sp.latex(F))" |> display
    display(mk.parse("---"))
    F = F(f(S)=>log(S)).doit().simplify()
    F |> display
    #FS = FS.removeO()(dS=>(μ*dt + σ*dW)*S).doit()    
    #display(FS)
end

$f{\left(S \right)} + ΔS \left. \frac{d}{d \xi_{1}} f{\left(\xi_{1} \right)} \right|_{\substack{ \xi_{1}=S }} + \frac{ΔS^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}} f{\left(\xi_{1} \right)} \right|_{\substack{ \xi_{1}=S }}}{2} + O\left(ΔS^{3}\right)$

$\frac{ΔS^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}} f{\left(\xi_{1} \right)} \right|_{\substack{ \xi_{1}=S }}}{2} + ΔS \left. \frac{d}{d \xi_{1}} f{\left(\xi_{1} \right)} \right|_{\substack{ \xi_{1}=S }} + f{\left(S \right)}$

$\frac{ΔS^{2} \frac{d^{2}}{d S^{2}} f{\left(S \right)}}{2} + ΔS \frac{d}{d S} f{\left(S \right)} + f{\left(S \right)}$

$\frac{S^{2} \left(ΔW σ + Δt μ\right)^{2} \frac{d^{2}}{d S^{2}} f{\left(S \right)}}{2} + S \left(ΔW σ + Δt μ\right) \frac{d}{d S} f{\left(S \right)} + f{\left(S \right)}$

$\frac{S^{2} \left(ΔW σ + Δt μ\right)^{2} \frac{d^{2}}{d S^{2}} f{\left(S \right)}}{2} + S \left(ΔW σ + Δt μ\right) \frac{d}{d S} f{\left(S \right)}$

$f(S+ΔS)-f(S)=\frac{S^{2} \left(ΔW σ + Δt μ\right)^{2} \frac{d^{2}}{d S^{2}} f{\left(S \right)}}{2} + S \left(ΔW σ + Δt μ\right) \frac{d}{d S} f{\left(S \right)}$

$ΔW σ + Δt μ – \frac{\left(ΔW σ + Δt μ\right)^{2}}{2}$

Ito Lemma¶

$$ \frac{dS_t}{S_t} = \mu dt + \sigma dW_t\\ dS_t = \left(\mu dt + \sigma dW_t\right)S_t $$

In [9]:

sp.@vars x X μ Z dt dW t W S dS S₀ σ T
sp.@symfuns f

F = f(S+dS)
FS = f(S+dS).series(dS,n=3)
L"f(S+dS)=%$(sp.latex(FS))" |> display
FS = FS.removeO().doit()
L"f(S+dS)=%$(sp.latex(FS))" |> display
L"dS=(\mu dt + \sigma dW)S" |> display
# 반드시 doit()  그래야 계수(편미분)이 S로 치환이 됨
FS = FS(dS=>(μ*dt + σ*dW)*S).simplify()
L"f(S+dS)=%$(sp.latex(FS))" |> display
   
dFS = FS - f(S)
dFS = dFS.doit()
L"df(S)=%$(sp.latex(dFS))" |> display

dFS = dFS.expand().doit()
# (dW^2=>dt,dW=>0,dW*dt=>0,dt^2=>0).
L"df(S)=%$(sp.latex(dFS))" |> display
L"(dW)^2=dt, dW\cdot dt=0, (dt)^2 = 0" |> display
dFS = dFS(dW^2=>dt,dW*dt=>0,dt^2=>0).doit().collect(dt)
L"df(S)=%$(sp.latex(dFS))" |> display

$f(S+dS)=f{\left(S \right)} + dS \left. \frac{d}{d \xi_{1}} f{\left(\xi_{1} \right)} \right|_{\substack{ \xi_{1}=S }} + \frac{dS^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}} f{\left(\xi_{1} \right)} \right|_{\substack{ \xi_{1}=S }}}{2} + O\left(dS^{3}\right)$

$f(S+dS)=\frac{dS^{2} \frac{d^{2}}{d S^{2}} f{\left(S \right)}}{2} + dS \frac{d}{d S} f{\left(S \right)} + f{\left(S \right)}$

$dS=(\mu dt + \sigma dW)S$

$f(S+dS)=\frac{S^{2} \left(dW σ + dt μ\right)^{2} \frac{d^{2}}{d S^{2}} f{\left(S \right)}}{2} + S \left(dW σ + dt μ\right) \frac{d}{d S} f{\left(S \right)} + f{\left(S \right)}$

$df(S)=\frac{S^{2} \left(dW σ + dt μ\right)^{2} \frac{d^{2}}{d S^{2}} f{\left(S \right)}}{2} + S \left(dW σ + dt μ\right) \frac{d}{d S} f{\left(S \right)}$

$df(S)=\frac{S^{2} dW^{2} σ^{2} \frac{d^{2}}{d S^{2}} f{\left(S \right)}}{2} + S^{2} dW dt μ σ \frac{d^{2}}{d S^{2}} f{\left(S \right)} + \frac{S^{2} dt^{2} μ^{2} \frac{d^{2}}{d S^{2}} f{\left(S \right)}}{2} + S dW σ \frac{d}{d S} f{\left(S \right)} + S dt μ \frac{d}{d S} f{\left(S \right)}$

$(dW)^2=dt, dW\cdot dt=0, (dt)^2 = 0$

$df(S)=S dW σ \frac{d}{d S} f{\left(S \right)} + dt \left(\frac{S^{2} σ^{2} \frac{d^{2}}{d S^{2}} f{\left(S \right)}}{2} + S μ \frac{d}{d S} f{\left(S \right)}\right)$

GBM (Generic Brownian Motion) Model¶

https://sine-qua-none.tistory.com/36

In [10]:

L"df(S)=%$(sp.latex(dFS))" |> display
L"f(S)=log(S)" |> display
df = dFS(f(S)=>log(S)).doit().collect(dt)
L"dlog(S)=%$(sp.latex(df))" |> display

$df(S)=S dW σ \frac{d}{d S} f{\left(S \right)} + dt \left(\frac{S^{2} σ^{2} \frac{d^{2}}{d S^{2}} f{\left(S \right)}}{2} + S μ \frac{d}{d S} f{\left(S \right)}\right)$

$f(S)=log(S)$

$dlog(S)=dW σ + dt \left(μ – \frac{σ^{2}}{2}\right)$

case1. $t=0$과 미래시점 $t=T$사이의 적분¶

여기 결과를 아래 $S_T$를 계산하는 데 사용

In [11]:

sp.@vars s w
S_T = sp.symbols("S_T")
# W₀ = 0
# T > 0
L"\int_{S_0}^S{log(s)}dz = \int_{W_0=0}^W{log(z)}dz" |> display
F = sp.integrate(1/s,(s, S₀, S)) - sp.integrate(μ-σ^2/2,(t,0,T)) - sp.integrate(σ,(w,0,W))
F |> display
S_T = sp.solve(F,S)[1]    
L"S_T=%$(latexify(S_T,env=:raw))" |> display

$\int_{S_0}^S{log(s)}dz = \int_{W_0=0}^W{log(z)}dz$

$- T \left(μ – \frac{σ^{2}}{2}\right) – W σ + \log{\left(S \right)} – \log{\left(S₀ \right)}$

$S_T=S_0 \cdot e^{T \cdot \mu – \frac{T \cdot \sigma^{2}}{2} + W \cdot \sigma}$

case2. 어떤(짧은 인터벌) $t=t_1$ 과 미래 시점 $t=t_2$ 사이의 적분¶

In [12]:

let
    sp.@vars T S₀ t₁ t₂ Wₜ₁ Wₜ₂ Sₜ₁ Sₜ₂
    sp.@vars S t W
    F = sp.integrate(1/S,(S, Sₜ₁,Sₜ₂)) - sp.integrate(μ-σ^2/2,(t,t₁,t₂)) - sp.integrate(σ,(W,Wₜ₁,Wₜ₂))
#     F = F.collect(t₂- t₁,sp.factor,exact=true)
    F = sp.sympy.powsimp(F,force=true)
    display(F)
    SOL = sp.solve(F,Sₜ₂)[1]
    display(sp.sympy.powsimp(SOL))
    SOL = SOL(t₂=> t₁+ 1//365).doit().simplify()
    
end

$Wₜ₁ σ – Wₜ₂ σ + t₁ \left(μ – \frac{σ^{2}}{2}\right) – t₂ \left(μ – \frac{σ^{2}}{2}\right) – \log{\left(Sₜ₁ \right)} + \log{\left(Sₜ₂ \right)}$

$Sₜ₁ e^{- Wₜ₁ σ + Wₜ₂ σ – t₁ μ + \frac{t₁ σ^{2}}{2} + t₂ μ – \frac{t₂ σ^{2}}{2}}$

Out[12]:

$Sₜ₁ e^{- Wₜ₁ σ + Wₜ₂ σ + \frac{μ}{365} – \frac{σ^{2}}{730}}$

sympy.stats : https://github.com/JuliaPy/SymPy.jl/issues/262

만기시점에 $S_T$ 생성¶

만기시점에 예상되는 만기 주가 $S_T$를 10,000개로 히스토그램 작성 및 만기 주가의 평균값이 얼마나 오르는지 확인
1번처럼 만기 주가를 10,000개 뽑으면서 만기 주가의 평균이 계속 달라지고 업데이트 되는 과정 관찰
$S_T$를 계산하는데 있어
- Case 1 : SymPy와 SymPy의 lambdify를 이용하는 경우
- Case 2 : 직접 코딩하는 경우
- Case 3 : SymPy만 이용하는 경우

$S_T$ 의 기대값 $\mathbb{E}\left(S_T\right)$, 분산 $\mathbb{V}\left(S_T\right)$, PDF 구하기¶

In [13]:

L"$S_T=$%$(latexify(S_T))" |> display
mk.parse("---") |> display

# W = √T ⋅ N(0,1)
ST = S_T(W=>√T*sp.Normal(x,0,1))

# Sₜ 의 기대치 E(Sₜ)
E = sp.E(ST)
theoretical_E = E.simplify()
L"$\mathbb{E}(S_T)=$ %$(latexify(theoretical_E))" |> display
mk.parse("---") |> display

# Sₜ의 분산 V(Sₜ)
V = sp.variance(ST)
theoretical_V = V.simplify()
L"$\mathbb{V}(S_T)=$ %$(latexify(theoretical_V))" |> display

# Sₜ의 표준편차 std(Sₜ)
σₜ = sp.std(ST)
theoretical_σₜ = σₜ.simplify()
L"$\mathbb{\sigma_T}(S_T)=$ %$(latexify(theoretical_σₜ))" |> display
mk.parse("---") |> display

# Sₜ 의 Probability Density Function (PDF)

# sp.density(St)의 결과 값은 Lambda 포맷으로 sp.Lambda를 사용하여 
# 결과를 사용할 수 있도록 한다.
theoretical_ST_PDF = sp.Lambda(sp.density(ST))(x)
L"\text{Probability Dendisty Function of}\; S_T \;:\; %$(latexify(theoretical_ST_PDF,env=:raw))" |> display

mk.parse("---") |> display
# PDF lambdify
# lamdify과정 : sp.Lambda(..)(x) =>  lambdify((..),(변수목록))
ST_PDF =  sp.lambdify( sp.Lambda(sp.density(ST))(x),(S₀,μ,σ,T, x))

z = 0:300
# StValues = convert.(Float64,StPDF.(100,0.02,0.2,1,z))
ST_Values(z) = ST_PDF.(100,0.02,0.2,1,z)    

L"S_T \text{ 는 Lognormal 분포를 따름-}\log{S_T}
        \text{가 정규분포를 따르기 때문에 }
        S_T\text{는 Lognormal 분포를 따른다.}" |> display

plot(z,ST_Values(z), title=L"Probability Density Function of $S_T$", titlefont=font(10),
    legend=false)

$S_T=$$S_0 \cdot e^{T \cdot \mu – \frac{T \cdot \sigma^{2}}{2} + W \cdot \sigma}$

$\mathbb{E}(S_T)=$ $S_0 \cdot e^{T \cdot \mu}$

$\mathbb{V}(S_T)=$ $S_0^{2} \cdot \left( e^{T \cdot \sigma^{2}} – 1 \right) \cdot e^{2 \cdot T \cdot \mu}$

$\mathbb{\sigma_T}(S_T)=$ $\sqrt{\left( – S_0^{2} \right) \cdot \left( 1 – e^{T \cdot \sigma^{2}} \right) \cdot e^{2 \cdot T \cdot \mu}}$

$\text{Probability Dendisty Function of}\; S_T \;:\; \frac{\sqrt{2} \cdot e^{\frac{ – \left( \left( – T \right) \cdot \mu + \frac{T \cdot \sigma^{2}}{2} + \log\left( \frac{x}{S_0} \right) \right)^{2}}{2 \cdot T \cdot \sigma^{2}}}}{2 \cdot \sqrt{pi} \cdot \mathrm{Abs}\left( \sqrt{T} \cdot x \cdot \sigma \right)}$

$S_T \text{ 는 Lognormal 분포를 따름-}\log{S_T} \text{가 정규분포를 따르기 때문에 } S_T\text{는 Lognormal 분포를 따른다.}$

Out[13]:

테스트 시나리오¶

조건이 아래와 같을 때 만기시점($T$)의 주가($S_T$)를 $n_{sample}$ 개 만큼 구하고 기대치 $\mathbb{E}\left(S_T\right)$ 계산

$S_0 = 100$ : 현재 시점 주가
$\mu=0.02$ : 주식의 연 수익률 $2\%$
$\sigma = 0.2$ : 변동성(주식의 연간 수익률의 표준편차)
$T=1$ : 만기시점 1년
$S_T$ : 만기시점의 주가

In [14]:

let    
    
    # sp.free_symbols(theoretical_E) |> display
    # sp.free_symbols(theoretical_σₜ) |> display
    E = sp.lambdify(theoretical_E,(S₀,μ,T))
    σₜ= sp.lambdify(theoretical_σₜ,(S₀,μ,σ,T))
    
    t=1:0.1:100
    s0,mu,sigma = 100, 0.02, 0.2
    
    
    p1 = plot(t,E.(s0,mu,t),
            title=L"\mathbb{E}\left(S_T\right)",
            xlabel=L"S_T",font=font(5),        
            titlefont=font(10),xguidefontsize=7
            )
    
    p2 = plot(t,σₜ.(s0,mu,sigma,t),
            title=L"{\sigma_T}\left(S_T\right)",
            xlabel=L"S_T",font=font(5),        
            titlefont=font(10),xguidefontsize=7
            )

    plot(p1,p2,size=(800,400),layout=2,legend=false)
    
end

Out[14]:

Case 1 : SymPy symbolic function 을 SymPy lambdify로 Julia function으로 변환 하여 사용¶

(아주 빠름)¶

장점 : symbolic function을 손으로 julia function으로 코딩 하지 않고 lambdify로 변환하기 아주 편리함

In [15]:

function make_GBM_path_1(;s0=100,mu=0.02,sigma=0.2,t=1,n_sample=1000,seed=nothing)
    
    !isnothing(seed) && Random.seed!(seed)
    
    z_normal = randn(n_sample)
    
    fn_st = sp.lambdify(S_T,(S₀,μ,σ,T,W))
    st = @. fn_st(s0,mu,sigma,t,√t * z_normal)
    
    # 이론 기대값
    fn_theoretical_avg = sp.lambdify(theoretical_E,(S₀,μ,T))
    theoretical_avg = fn_theoretical_avg(s0,mu,t)
    
            
    # 이론 PDF     
    ST_PDFValues(z) = ST_PDF.(s0,mu,sigma,t, z)

    visualize(st,theoretical_avg,n_sample,ST_PDFValues)
end
#  @time make_GBM_path_1(seed=123)

Out[15]:

make_GBM_path_1 (generic function with 1 method)

Case 2 : SymPy symbolic function 을 직접 코딩하여 Julia function으로 변환 하여 사용¶

(최고 빠름)¶

단점 : symbolic function을 손으로 julia function으로 코딩해야 하는 번거로움 있음

In [16]:

function make_GBM_path_2(;s0=100,mu=0.02,sigma=0.2,t=1,n_sample=1000,seed=nothing)
    
    !isnothing(seed) && Random.seed!(seed)
    
    z_normal = randn(n_sample)
    
    drift = (mu - 0.5 * sigma^2)*t
    diffusion = sigma * √t
    
    st = @. s0 * exp(drift + diffusion * z_normal)
    
    # 이론 기대값
    theoretical_avg = s0*exp(mu*t)
    
    # 이론 PDF     
    ST_PDFValues(z) = ST_PDF.(s0,mu,sigma,t, z)
    
    visualize(st,theoretical_avg,n_sample,ST_PDFValues)
end
# @time make_GBM_path_2(seed=123)

Out[16]:

make_GBM_path_2 (generic function with 1 method)

Case 3 : SymPy symbolic function 을 사용¶

(가장 느림)¶

In [17]:

function make_GBM_path_3(;s0=100,mu=0.02,sigma=0.2,t=1,n_sample=1000,seed=nothing)
    
    !isnothing(seed) && Random.seed!(seed)
    
    z_normal = randn(n_sample)
    
    # z_normal 이 vector로써 Sₜ를 vectorize 해야 함
    # Sₜ 리턴값은 symbolic 술자 목록 𝑆𝑦𝑚𝑃𝑦.𝑆𝑦𝑚[] 으로 N()으로
    # 숫자로 변환 한다.
    
    st = sp.N(@. S_T(S₀=>s0,μ=>mu, σ=>sigma, T=>t, W=> √t*z_normal))
    
    # 이론 기대값
    theoretical_avg = sp.N(theoretical_E(S₀=>s0, μ=>mu, T=>t))
    
    # 이론 PDF     
    ST_PDFValues(z) = ST_PDF.(s0,mu,sigma,t, z)
    
    visualize(st,theoretical_avg,n_sample,ST_PDFValues)
end
# @time make_GBM_path_3(seed=123)

Out[17]:

make_GBM_path_3 (generic function with 1 method)

Visualization¶

In [18]:

function visualize(st,theoretical_avg,n_sample,ST_PDFValues)
    # Sₜ 기대값
    st_mean = mean(st)
    st_step_avg = map(i->mean(st[1:i]),2:n_sample)    
    # Sₜ PDF
    z = 0:300
    st_pdf = ST_PDFValues(z)
    
    L"""
    W \sim \sqrt{T}\cdot N(0,1)\\
    \text{Average:}S_T\;
    \mathbb{E}\left(S_T\right)= 
    \mathbb{E}\left( %$(latexify(S_T,env=:raw))\right)=
    %$(round(st_mean,digits=4))\\
    
    \text{Theoretical average:}\;\mathbb{E}(S_T)=
    %$(latexify(theoretical_E,env=:raw))=
    %$(round(theoretical_avg,digits=4))
    """|> display
    
    mk.parse("---") |> display
    plot(size=(800,300),layout=2)
    stephist!(st, title=L"Chart of $S_T$",
        #legend=false,
        titlefont = font(10),normalize=true,
        fill=true, fillalpha=0.5,
        label=L"$S_T$ histogram",
        subplot=1)
    vline!([st_mean],color=:red,lw=3,
        label=L"$S_T$ mean value",
        subplot=1)
    plot!(z,st_pdf,label=L"$S_T$ PDF",lw=2, color=:green,
        subplot=1)
    plot!(st_step_avg,layout=2,st=:line,
        title=L"Convergence of Average of $S_T$",
        legend=false,titlefont = font(10),
        subplot=2)    
end

Out[18]:

visualize (generic function with 1 method)

3 Case 경우 성능 시험¶

In [24]:

seed = Int(floor(datetime2unix(now())))

Out[24]:

1673258865

In [25]:

L"\text{Case 1 : SymPy lamdify 사용 (아주빠름)}" |> display
@time make_GBM_path_1(n_sample=10_000,seed=seed)

$\text{Case 1 : SymPy lamdify 사용 (아주빠름)}$

$W \sim \sqrt{T}\cdot N(0,1)\\ \text{Average:}S_T\; \mathbb{E}\left(S_T\right)= \mathbb{E}\left( S_0 \cdot e^{T \cdot \mu – \frac{T \cdot \sigma^{2}}{2} + W \cdot \sigma}\right)= 102.1624\\ \text{Theoretical average:}\;\mathbb{E}(S_T)= S_0 \cdot e^{T \cdot \mu}= 102.0201 $

  0.228980 seconds (363.90 k allocations: 396.789 MiB, 17.75% gc time, 46.75% compilation time)

Out[25]: