Hatena::ブログ(Diary)

My Life as a Mock Quant このページをアンテナに追加 RSSフィード

2018-04-23

N個の独立な一様乱数(-1, 1)の和の標準偏差は1/3√N

| 20:37 | N個の独立な一様乱数(-1, 1)の和の標準偏差は1/3√N - My Life as a Mock Quant を含むブックマーク

手計算で確認&ボケ防止にRでチェック。

f:id:teramonagi:20180423203650p:image

sizes <- 2**(1:15)
y <- numeric(length(sizes))
for(i in seq_along(sizes)){
  x <- numeric(10**3)
  for(j in seq_len(10**3)){
    x[j] <- sum(runif(sizes[i], min=-1, max=1))
  }
  y[i] <- sd(x)
}
plot(sizes, y)
lines(sizes, sqrt(1/3*sizes))

2018-03-30

purrr::map_dfr = lapply + dplyr::bind_rows

| 22:54 | purrr::map_dfr = lapply + dplyr::bind_rows - My Life as a Mock Quant を含むブックマーク

そういうことなんだよな〜シミュレーション系でよく使うのでメモ。

> dplyr::bind_rows(lapply(1:3, function(x){head(iris, 1)}))
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          5.1         3.5          1.4         0.2  setosa
3          5.1         3.5          1.4         0.2  setosa
> purrr::map_dfr(1:3, ~ head(iris, 1))
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          5.1         3.5          1.4         0.2  setosa
3          5.1         3.5          1.4         0.2  setosa

2017-01-15

確率変数を変数変換した場合の確率分布

| 15:57 | 確率変数を変数変換した場合の確率分布 - My Life as a Mock Quant を含むブックマーク

昔やった気もするが、すぐ忘れる&頭の体操もかねてもう一度

算数

適当な確率分布に従う確率変数x(面倒なので[0, 1]区間の一様分布想定)とそれを2乗した変数yを考える。

このときyの従う確率分布は

 1 = ¥int_0^1 p(x) dx = ¥int_0^1 1 dx = ¥int_0^1 1/(2*¥sqrt{y}) dy

より

 p(y) = 1/(2¥sqrt{y})

となる。

Rでやる

0, 1区間の一様な乱数x, 及びそれを二乗したyを生成。

x <- runif(10^5, 0, 1)
y <- x^2

yのヒストグラムが↑の計算と合うかチェック

hist(y, probability=TRUE)
t <- seq(0, 1, by=0.01)
lines(t, 1/(2*sqrt(t)),col="red", lwd=3)

f:id:teramonagi:20170115155454p:image

2016-08-16

金融市場が混乱する間隔を日経平均ボラティリティー・インデックス(日経平均VI)から推定したい

| 18:48 | 金融市場が混乱する間隔を日経平均ボラティリティー・インデックス(日経平均VI)から推定したい - My Life as a Mock Quant を含むブックマーク

日経ボラティリティー・インデックス(日経平均VI)祭りはまだ続く。

やりたいことは簡単で

  • 日経VIが跳ねるタイミング=金融市場が混乱

だと考えて、"日経VIが跳ねるタイミング"の確率分布がどのような形になるのかを推定するというお話。

個人的には指数分布なんじゃね?と思っていたが、違った。

データの準備

ほぼ

でやった内容とおなじで、日経VIが跳ねたポイント(下図赤点)を手でやるのがめんどうなので、Rに同定させるためのコード。

↑のLINKとはrollapp()関数がちょっと違ってるので注意

library("Quandl")
# Sorce: https://www.quandl.com/data/NIKKEI/VLTL-Nikkei-Stock-Average-Volatility-Index
# Download Nikkei VI
vi <- Quandl("NIKKEI/VLTL")
#Decay detect
rollapp <- function(x, k, fun)
{
  w <- k/2
  sapply(seq_len(length(x)), function(i){
    index_s <- max(i-w, 1)
    index_e <- min(i+w, length(x))
    middle <- index_s + i-w
    fun(1+i-index_s, x[index_s:index_e])
  })
}
x <- vi$Close
max_point <- rollapp(x, 74, function(i, x){x[i]==max(x)})
max_point <- replace(max_point, is.na(max_point), FALSE)
plot(x, type="l")
points(seq_len(length(x))[max_point], x[max_point], col = "red", pch=16)

f:id:teramonagi:20160816185327p:image

このデータの赤点同士の時間間隔がどういう確率分布に従っているのかを評価したい。

推定(指数分布)

指数分布で確率密度関数を推定してみる

# 適当にdata.frame化
df <- dplyr::group_by(data.frame(index=cumsum(max_point)), index) %>% summarize(n=n())
lambda <- 1/mean(df$n)
# 下記でも同じ
# library("MASS")
# exponential <- coef(fitdistr(df$n, "exponential"))
hist(df$n, freq=FALSE)
curve(dexp(x, rate=lambda),add=TRUE,lty=2)

f:id:teramonagi:20160816185326p:image

う、うーん、これではない感。

推定(ワイブル分布)

によると、ワイブル分布は「時間に対する劣化現象や寿命を統計的に記述するためにも利用される」とあるので、こいつでやってみたい。

金融市場が破裂する=寿命が尽きると考えるわけだ。

library("MASS")
weibull <- coef(fitdistr(df$n, "weibull"))
hist(df$n, freq=FALSE)
curve(dweibull(x, weibull["shape"], weibull["scale"]), add=TRUE, lty=2)

f:id:teramonagi:20160816185325p:image

これだなー!

あとはこの裏側の寿命強度過程的なものがどうなっているのかを調べれば大体OKだな。

参考

2015-11-08

確率変数の変換について

| 21:24 | 確率変数の変換について - My Life as a Mock Quant を含むブックマーク

ロジックメモ

ある確率変数 X: ¥Omega ¥rightarrow ¥mathbb{R}とその写像 Y=g(X): ¥Omega ¥rightarrow ¥mathbb{R}で定義される2つの確率変数、特に確率変数 Yの確率分布について考えたい。

 Yの確率分布関数 F_Y(y)

 F_Y(y) = P(Y ¥le y)= P(g(X) ¥le y) =  P(¥left¥{ ¥omega ¥in ¥Omega; g(X(¥omega)) ¥le y¥right¥}) = P(¥left¥{ ¥omega ¥in ¥Omega; X(¥omega) ¥le g^{-1}(y)¥right¥}) = P(X ¥le g^{-1}(y)) = F_X(g^{-1}(y))

と確率変数 Xの確率分布関数を用いてあらわすことができる。ここで、 gは単調増加な関数であると仮定している。単調減少の場合は不等号の向きが逆になる。それ以外の場合は、単調増加・現象とみなせる区間に切って、確率を分けて計算し、最後に足しあげる操作が必要。

したがって、確率変数 Yの確率密度関数f_Yは確率変数 Xの確率密度関数f_Xを用いて

 f_Y(y) = ¥frac{d}{dy} F_Y(y) = f_X(g^{-1}(y)) ¥cdot ¥frac{d}{dy} g^{-1}(y)

と書くことができる。

Rで確かめる

これを確かめてみよう。ここでは g(x) = ¥log x^2と定義する。すると上式から

 F_Y(y) = P(¥left¥{ ¥omega ¥in ¥Omega; ¥log X(¥omega)^2 ¥le y¥right¥}) = P(¥left¥{ ¥omega ¥in ¥Omega; - ¥exp(y/2) ¥le X(¥omega) ¥le ¥exp(y/2) ¥right¥}) = F_X(¥exp(y/2)) - F_X(-¥exp(y/2))

とすることができる。ここで関数 gx=0を境に単調減少から単調増加関数へと変わるので、上のように項が2つ出てくる。上で言及した最後の"足しあげ"に相当するものだ。

従って確率密度関数 f_Y(y)

 f_Y(y) = f_X(¥exp(y/2)) ¥cdot ¥frac{1}{2}¥exp(y/2) - f_X(-¥exp(y/2)) ¥cdot ¥frac{1}{2}(-¥exp(y/2)) = ¥frac{1}{2}¥exp(y/2) ¥left¥{ f_X(¥exp(y/2)) + f_X(-¥exp(y/2))¥right¥}

となる。更に(本来、何でもいいんだが)確率変数 Xが標準正規分布に従うと仮定して上の結果を確かめてみる。

size <- 10^5
#標準席分布に従う乱数の生成
x <- rnorm(size)
#関数gで変換
y <- log(x^2)

この結果をプロットする(確率密度表示)&上で計算した確率密度関数を上から重ねてみると

hist(y, sqrt(size), freq=FALSE)
ys <- seq(-20, 5, 0.1)
fy <- function(y){exp(y/2)/2*(dnorm(exp(y/2)) + dnorm(-exp(y/2)))}
lines(ys, fy(ys), col=2, lwd=3)

f:id:teramonagi:20151108212119p:image

・・・となってほぼぴったり一致する。なので、計算は正しいだろうと。

ggplot2に慣れたいんで、ggplot2でもやってみた。..density..をggplot関数じゃなくてgeom_histogram関数の中に書かなきゃだめなのがミソなんだが、こういうの覚えきれん。。。

library("ggplot2")
ggplot(data.frame(value=y), aes(x=value)) + 
  geom_histogram(aes(y=..density..), binwidth=0.1, fill="light blue") + 
  stat_function(fun=fy, colour="blue", size=2)

f:id:teramonagi:20151108221617p:image