RでOR:金融工学

目次

  1. オプション価格
    1. ブラックショールズ公式
    2. インプライドボラティリティ
    3. オプション価格の2項モデル公式
  2. 確率微分方程式
    1. ブラウン運動
    2. t分布でランダムウォーク
    3. 幾何ブラウン運動
    4. オプションの価格、非ブラウン運動の場合
  3. 金利モデル
    1. Vasicekモデル
    2. Cox-Ingasol-Rubinsteinモデル
  4. データ分析
    1. 日経225平均株価の収益率の自己相関
    2. 収益率の絶対値(平均絶対偏差)の分析
  5. ARCH型モデル
    1. ARCHモデルのシミュレーション
    2. GARCHモデルのシミュレーション
    3. ボラティリティクラスタリング検証実験

リスク資産価値の確率変動モデル

次の図は東京株式市場の日経平均株価終値の推移を時系列グラフで表したものである。このようなランダムな動きをモデル化したものにランダムウォークやブラウン運動がある。

    n225

ランダムウォーク

2種類の数 +1, -1 をランダムに生成し、その部分和のグラフを描くと、左右にランダムに歩いた軌跡のようなものができることから、これをランダムウォークと呼ぶ。すなわち、2値 +1, -1 を等確率で取る独立な確率変数列 {Xn} の部分和を {Sn}としたとき、{Sn} をランダムウォークという。

Sn = Sn-1 + Xn, n = 1,2,...; S0 = 0
P(Xn = 1) = P(Xn = -1) = 0.5, n=1,2,...

Sn は、パラメータ n, p=0.5 の2項分布に従うので、その変動幅は n の平方根のオーダーである。次のプログラムはそのようなランダムウォークのサンプルパスを重ね描きして、±2√n のグラフを付加するものである。実行結果例を添える。

> n = 400
> plot(c(0,n), c(-2.5,2.5)*sqrt(n), type="n") > for(i in 1:10) lines(c(0, cumsum(sample(c(-1,1), n, replace=TRUE))), col=i) > abline(h=0) > curve(-2*sqrt(x), add=TRUE, lty=2) > curve(2*sqrt(x), add=TRUE, lty=2)

    random

ランダムウォークを微小時間 Δt 間隔で繰り返すことを考える。Xn の取りうる値も小さくしないとおかしな動きになるので ±Δa とする。区間 [0,t] の間に n 回推移したとして(t = nΔt)t におけるランダムウォークの値を S(t) とすると、S(t) はドモアブルラプラスの定理により、平均0分散 nΔa2 = tΔa2 /Δt の正規分布に従う。Δt をゼロに近づけても分散が有限の範囲で留まるために、Δa = σ√(Δt) を仮定すると、結局、S(t) は近似的に平均0分散 σ2t の正規分布に従う。Δt →0 とした時の極限の確率過程を B(t) と書いてブラウン運動と呼ぶ。

ランダムウォークのグラフと同じような例を載せる。Δt=0.001 とそれほど小さくないが、軸の目盛りを無視すれば、2つの違いはほとんど感じられれない。

> t = 4
> n = 100
> plot(c(0,t), c(-2.5,2.5)*sqrt(t), type="n")
> for(i in 1:10) lines(seq(0,t,1/n), c(0, cumsum(rnorm(n*t,0,sqrt(1/n)))), col=i)
> abline(h=0)
> curve(-2*sqrt(x), add=TRUE, lty=2)
> curve(2*sqrt(x), add=TRUE, lty=2)

    brown

リスク資産価格に基づく数量分析ではブラウン運動を仮定することが多い。

負の値を取らないようにするために、その対数を取ったものがブラウン運動をするという確率過程もよく使われる。これは幾何ブラウン運動と呼ばれる。

オプション価格

リスク商品のコールオプションは、将来のある時点でその商品をある価格で購入することができる権利のことである。プットオプションは、売ることができる権利のことである。コールオプション、プットオプションはもとのリスク商品の派生証券と呼ばれ、もとのリスク商品は原証券と呼ばれる。

権利行使時点を満期日、権利行使時に決済する価格を(権利)行使価格あるいはストライクプライス、という。時点 t における原証券の価格を S(t)、満期日を T、権利行使価格を K とすると、コールオプションの価値は max{S(T) - K, 0} と表される。従って、その派生証券の現在価値(t = 0 での価格、S(0) は既知)は、安全利子率を r とすると、e-rT E(max{S(T) - K, 0} となる。S(T) の分布がわかれば、現在価値を計算することができる。

例えば、S(T) が平均 S(0)erT、分散 Ts2 の正規分布に従っているとすれば、正規乱数を使って、現在価格を数値的に求めることができる。

ブラックショールズ価格式

ブラックショールズ価格式(あるいは、公式、Black-Scholes formula)とは、原証券の価格過程が幾何ブラウン運動に従うと仮定して、与えられた初期価格(S)、ボラティリティ(s)、満期日(T)、行使価格(K)、安全利子率(r) に対して、コールオプションの価格を求めるものである。

関数化したプログラムを示す。ここで、pnorm 関数は標準正規分布の累積確率を返す。

blackScholes = function(K, s, T, r, S=100) {
	d1 = (log(S/K) + (r+(s^2)/2) * T) / s / sqrt(T)
	d2 = (log(S/K) + (r-(s^2)/2) * T) / s / sqrt(T)
	return(S * pnorm(d1) - K * exp(-r*T) * pnorm(d2))
}

# 実行例:満期の違いを比較(T の単位は年)
> r = 0.03; K = 110; T = 0.25; s = 0.2; S = 100 > blackScholes(K, s, T, r, S) [1] 1.091344 > blackScholes(K, s, c(0.25,0.5,1), r, S) [1] 1.091344 2.611902 5.293398

行使価格の違いによるオプション価格の変化を調べる。満期が直後に迫っていれば、オプション価格はほぼ max(0, S-K) に等しいと考えられる。満期日が先になるにつれて、債権価格変動が大きくなり、オプションの時間価値が増えるので、このラインから上にシフトする。これらの定性的な性質を数量的に検証する。

次がそのプログラムの例である。

> strike = 80:120
> T = 0.25
> bs = blackScholes(strike, s, T, r, S)
> curve(ifelse(S > x, S-x, 0), xlim=c(80,120), lty=2, asp=1) > abline(h=0) > lines(strike, bs, type="l")
> T = 0.5 > bsa = blackScholes(strike, s, T, r, S) > lines(strike, bsa, type="l", col=3) > legend(100, 20, c("T=0.25", "T=0.5"), lty=1, col=c(1,3))

上で定義した blacksholes 関数は引数をベクトルとしても構わないので、計算したい行使価格のベクトルを行使価格変数として与えれば、結果はそれぞれの行使価格に対するオプション価格のベクトルになる。
curve(ifelse(... は、参考のために max(S-x,0) のグラフを描いたものである。curve 関数の最初の引数は関数名を指定するのが通常だが、このように関数定義を書いても構わない。ifelse(C, a, b) は条件 C が TRUE ならば a、さもなければ b を返す関数である。asp=1 は、漸近線の傾きを -45度に描くためのオプション指定である。
legend 関数は凡例を描く。

次は実行結果のグラフである。満期が延びるにつれて、不確実性の時間価値が増大し、オプション価格が上昇している様子がわかる。

    bs 


インプライドボラティリティ

ブラックショールズ公式に含まれるパラメータのうち、ボラティリティを除く4つ、初期価格、満期日、行使価格、安全利子率は確定しているが、ボラティリティについては、いろいろな推定法があり、一意に決まらない。ブラックショールズ公式に基づくオプション価格(BS価格と言おう)が与えられた時、その価格を導くために使われたボラティリティの値をインプライドボラティリティという。

インプライドボラティリティは、ボラティリティを未知の変数としてオプション価格が BS価格に等しいという方程式を、ボラティリティについて解いたもので、これは、方程式のゼロ点を求める問題となる。

ブラックショールズ公式をボラティリティの関数としてグラフを描いてみる。

> curve(blackScholes(K,x,T,r,S), 0, 0.4, bty="n")

curve 関数は x を独立変数とする関数のグラフを描くので、blackScholes 関数のように独立変数(パラメータ)がたくさんある場合は、グラフを描こうする変数を x として指定する必要がある。結果は次のようになる。

    impv

BS 価格が与えられた時のインプライドボラティリティは、上の関数グラフで BS価格を y切片とする水平線との交点で与えられる。関数の形状が単純なので、どんな方法でも容易に解を求めることができる。

Rの uniroot 関数は、連続関数で関数値の符号が異なる2点が見つかれば、その2点を端点とする閉区間に含まれるその関数のゼロ点を返す。引数の最初に関数名を書くが、多変数関数の場合は最初の変数に関するゼロ点を求めるという約束なので、ボラティリティを最初の変数にした関数を定義する必要がある。

> fbs = function(s,K,T,r,S,p) blackScholes(K,s,T,r,S) - p
> uniroot(fbs, c(0,1), p=5, K=110, T=0.5, r=0.03, S=100)
$root
[1] 0.2911707

$f.root
[1] -2.425306e-06

$iter
[1] 5

$estim.prec
[1] 8.363982e-05

オプション価格の2項モデル公式

binomialModel = function(K,s,n,u,d,r) {
	riskn = (r - d) / (u - d)
	C = 0
	for(i in 0:n) 
		C = C + dbinom(i,n,riskn) * max((1+u)^i * (1+d)^(n-i) * s - K, 0)
	return (C / (1+r)^n)
}
u = exp(0.2/sqrt(250))-1; d = 1/(1+u)-1; r = 0.05/250; s = 100; K = 105
n=30
z = sapply(c(90:110), binomialModel,s,n,u,d,r)
plot(c(90:110),z, type="l")
lines(c(s,90),c(0,s-90),col=2)

2項格子モデルのシ ミュレーション(CRR)

####### オプション評価の格子モデルシミュレーション ###################
n  = 12          # 1年の期間数
dt = 1/n         # 1期間の長さ(年単位)
m  = 5           # 満期までの期間
r  = 0.1         # 無リスク金利
a  = exp(r * dt) # 割引関数 
sd = 0.4         # 標準偏差
u  = exp(sd * sqrt(dt)); d = 1/u # 格子モデルの推移確率
p  = (a - d) / (u - d)           # リスク中立確率
S0 = 50          # 初期株価
K  = 50          # 行使価格
cat(n, dt, m, r, a, sd, u, d, p, S0, K)

#### 価格変化
x = matrix(rep(0,(m+1)*(m+1)), m+1)
x[1,1] = S0
for(i in 1:m) {
  for(j in 1:(i+1)) {
    x[i+1,j] = x[i,j] * d
  }
  x[i+1,i+1] = x[i,i] * u 
}
fix(x)

#### コールオプション価格
y = matrix(rep(0,(m+1)*(m+1)), m+1)
for(j in 1:(m+1)) {
  y[m+1,j] = max(0,x[m+1,j] - K)
}
for(i in m:1) {
  for(j in 1:i) {
    y[i,j] = max(((1-p) * y[i+1,j] + p * y[i+1,j+1]) / a, x[i,j] - K)
  }
}
fix(y)

#### プットオプション価格
z = matrix(rep(0,(m+1)*(m+1)), m+1)
for(j in 1:(m+1)) {
  z[m+1,j] = max(0,K - x[m+1,j])
}
for(i in m:1) {
  for(j in 1:i) {
    z[i,j] = max(((1-p) * z[i+1,j] + p * z[i+1,j+1]) / a, K - x[i,j])
  }
}
fix(z)

#### 一括出力(証券価格/コールオプション価格/プットオプション価格
w = matrix(rep(0, 3*(m+1)*(m+1)), 3*(m+1))
for(i in 1:(m+1)) {
  for(j in 1:i) {
    w[3*i-2,j] = x[i,j]
    w[3*i-1,j] = y[i,j]
    w[3*i-0,j] = z[i,j]
  }
}
fix(w)

ページトップにもどる 索引にもどる

確率微分方程式

ブラウン運動

################### ランダムウォーク ########################
m = 0; s = 0.1; dt = 0.1; n = 200
x = rnorm(n, m*dt, s*sqrt(dt)); x[1] = 1
plot(dt*c(0:(n-1)), cumsum(x), type="l", ylim=c(0.5,1.5))

t分布でランダム ウォーク

時々、大きな変位が観測される

m = 0; s = 0.1; dt = 0.1; df = 3; n = 200
x = rt(n,df) * s*sqrt(dt) + m*dt; x[1] = 1
plot(dt*c(0:(n-1)), cumsum(x), type="l", ylim=c(0.2,1.8), main=paste("df=",df))
x = rt(n,df) * s*sqrt(dt) + m*dt; x[1] = 1
lines(dt*c(0:(n-1)), cumsum(x))			# 重ねて描く

#### 比較するために、正規分布変位のランダムウォークを重ねて描く
y = rnorm(n, m*dt, s*sqrt(dt)); y[1] = 1
lines(dt*c(0:(n-1)), cumsum(y), col=2)
y = rnorm(n, m*dt, s*sqrt(dt)); y[1] = 1
lines(dt*c(0:(n-1)), cumsum(y), col=2)	

幾何ブラウン運動

対数を取ったものがランダムウォークするので、正規乱数の累積和(ランダムウォーク)を指数変換すればよい。

n = 200; m = 0; s = 0.01;
z = c(0,rnorm(n,m,s))
ss = cumsum(z)
plot(0:n, exp(ss), ylim=c(0.8,1.2), type="l")

### 一行で
lines(0:n, c(1,exp(cumsum(rnorm(n,m,s)))))

オプションの価格、 非ブラウン運動の場合

収益率がt分布に従う場合のオプション価格:収益率のパスを発生させて、max(S-K,0)の期待値を計算する

heavytail = function(r,s,df,S,T,dt) {
  n = T / dt
  za = cumsum(log(1 + dt*r+sqrt(dt)*s*rt(n,df)))   ### これがt乱数
  x = S * exp(za)
  ### plot(x, type="l")
  return(x[n])
}
df = 3; dt = 0.001; s = 0.2; r = 0.05; S = 100; K =105; T = 1
m = 1000; sec = sapply(1:m, function(i) heavytail(r,s,df,S,T,dt))
hist(sec)
for(i in 1:m) aaa[i] = max(sec[i]-K,0)
mean(aaa)

「rt」の代わりに「rnorm」を使えば、ブラックショールズ式のシミュレーションになる
かなり割高な価格付けになる。
自由度を小さくすると、logの中身がマイナスになる場合がある

もどる

金利モデル

Vasicekモ デル

初期値に回帰するように、ランダムウォークに初期値との差に比例する項を付け加えたものである。

vasicek = function(a, s, r0, n) {
  noise = rnorm(n)
  r = numeric(n)
  r[1] = r0
  for(t in 2:n) 
    r[t] = r[t-1] + a * (r0-r[t-1]) + s * noise[t]
  return(r)
}

# 実行例
> r = vasicek(0.05, 0.005, 0.05, 500) > plot(r, ylim=c(0,0.1), type="l", main=paste("a=",a,"s=",s)) > abline(h=r0, lty=2) > for(i in 2:5) { + r = vasicek(0.05, 0.005, 0.05, 500) + lines(r, col=i) + }

Cox-Ingasol -Rubinsteinモデル

Vasicek モデルのノイズ項の分散を状態に比例させて、ゼロに近づいた時に変動幅を抑える(金利がマイナスになることを防ぐ)工夫をしたものである。離散近似のため、平方根の中身を0で止めてある

CIR = function(a, s, r0, n) {
	r = numeric(n)
	r[1] = r0
	noise = rnorm(n)
	for(t in 2:n) 
		r[t] = r[t-1] + a * (r0 - r[t-1]) + s * sqrt(max(0,r[t-1])) * noise[t]
	return(r)  
}

# 実行例
> a = 0.01 > s = 0.01 > r0 = 0.05 > n = 500 > r = CIR(a, s, r0, n) > plot(r, ylim=c(0,0.10), type="l", main=paste("a=",a,"s=",s)) > abline(h=r0, lty=2) > for(i in 2:5) { + r = CIR(a, s, r0, n) + lines(r, col=i) + }

5通りのサンプルパスを描いた例を載せる。

    cir

もどる

データ分析

日経225平均株 価の収益率の自己相関

  1. 収益率は無相関
  2. 収益率の絶対値は強い正の相関
  3. 収益率の絶対値を収益率の標準偏差の推定値と考えると、標準偏差が一定ではない
n255 = read.table("clipboard")[,1] # Excelからデータをコピー
n = length(n255)

ret = abs(n255[2:n] - n255[1:(n-1)]) / n255[1:(n-1)] # 収益率の絶対値
n1 = n-1
r = sapply(1:50, 
     function(i) (mean(ret[1:(n1-i)] * ret[(i+1):n1]) - mean(ret)^2) / var(ret))
plot(r, type="l", ylim=c(-0.3,0.3)) 

ret = (n255[2:n] - n255[1:(n-1)]) / n255[1:(n-1)] # 収益率
n1 = n-1
r = sapply(1:50, 
     function(i) (mean(ret[1:(n1-i)] * ret[(i+1):n1]) - mean(ret)^2) / var(ret))
par(new=T)
plot(r, type="l", ylim=c(-0.3,0.3))               # 重ね描き

収益率の絶対値 (平均絶対偏差)の分析

日経平均終値の収益率を求め、その絶対値の移動平均を計算する

df = read.csv("stocktestdata.csv", header=T, skip=5)
final = df$終値; n = length(final)          # 終値だけ取り出す
rate = (final[-1] - final[-n]) / final[-n]  # 収益率
par(mfrow=c(2,1))
ts.plot(rate)            ## 収益率
ts.plot(abs(rate))       ## 収益率の絶対値

## 移動平均(sapply利用)
N = length(rate)
m = 10
MA.m = sapply(1:(N-m+1), function(i) mean(abs(rate[i:(i+m-1)])))
plot(abs(rate), type="l", col=3)
lines(c(rep(0,m-1),MA.m))

ページトップにもどる 索引にもどる

ARCH型モデル

ARCHモデルのシミュレーション

### ARCH モデルのシミュレーション

a0 = 0.5; a1 = 1 - a0;
n = 1000; x = c(0); sg = c(a0);
for(i in 1:n) {
	sg = c(sg, a0+a1*x[i]^2)
	x = c(x, rnorm(1)*sqrt(sg[i+1]))
}
plot(x, type="l", main=paste("ARCH (a0 a1=",a0,a1,")"))
plot(sqrt(sg), type="l", main="ボラティリティ")

win.graph()
par(mfrow=c(2,1))

acf(x^2)

df = read.csv("G:/stockdata.csv", header=T, skip=7); n = dim(df)[1]
rr = (df[-1,5]-df[-n,5])/df[-n,5]; m = length(rr)
plot(rr[4000:5000], type="l")
acf(abs(rr[4000:5000]))

GARCHモデルのシミュレーション

a0 = 0.5; a1 = 0.1; b = 0.85;
n = 1000; x = c(1); sg = c(a0);
for(i in 1:n) {
	sg = c(sg, a0+a1*x[i]^2+b*sg[i])
	x = c(x, rnorm(1)*sqrt(sg[i+1]))
}
plot(x, type="l", main=paste("GARCH (a0 a1 b=",a0,a1,b,")"))
plot(sqrt(sg), type="l", main="ボラティリティ")

ボラティリティクラスタリング検証実験

### 日経平均終値の収益率のボラティリティクラスタリング検証実験
### 収益率が独立同分布しているならば、異常値(絶対値がある値以上になる)の出現はポワソン過程に従うはず
### 異常値のギャップを求め、その間隔が指数分布に従っているかどうかをQ-Qプロットでチェックする

df = read.csv("G:/stockdata.csv", header=T, skip=7); n = dim(df)[1]
rr = (df[-1,5]-df[-n,5])/df[-n,5]; 
m = length(rr)
sr = sd(rr)

size = 3; gap = which(rr > size*sr); g = gap[-1] - gap[-length(gap)]
mg = mean(g); n = length(gap)-1

plot((c(1:n)-0.5)/n, pexp(sort(g),1/mg))
plot(sort(g), qexp((c(1:n)-0.5)/n,1/mg))

ページトップにもどる 索引にもどる