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

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

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))

purrr::map_dfr = lapply + dplyr::bind_rows

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

> 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

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

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

算数

適当な確率分布に従う確率変数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)

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

日経ボラティリティー・インデックス(日経平均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)

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

推定(指数分布)

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

# 適当に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)

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

推定(ワイブル分布)

によると、ワイブル分布は「時間に対する劣化現象や寿命を統計的に記述するためにも利用される」とあるので、こいつでやってみたい。
金融市場が破裂する=寿命が尽きると考えるわけだ。

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

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

確率変数の変換について

ロジックメモ

ある確率変数 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)

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

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)

売上げを気温で回帰してみる

あたまだし

掲題の件、最近出たわかりやすいデータ分析の入門書、

において

  • 売上げ v.s. 気温

として回帰分析をしていたケースがあった。これが時系列データに対するやってはいけない回帰、すなわち、みせかけの回帰に相当するケースではないか?という話を久保先生としていた。


ここでは、これを簡単なモデルのシミュレーションを通じてチェックしてみたい。

売上げ(平均値が正の正規分布)を気温(平均回帰過程)で回帰してみる

簡単のため、

としてモデリングする。コードはかつてないほどの殴り書き&モデルパラメータは適当。
kに対するforループが上記の過程に従う1パスのモンテカルロシミュレーションに対応し、それを何回も繰り返して、回帰して、気温に対する回帰係数のp値の分布を見てみる。

size_step <- 100
sigma <- 5
dt <- 0.1
theta <- 1
mu <- 25
p_value <- numeric(10^3)
for(k in 1:10^3)
{
  temp <- rep(mu, size_step+1)
  for(i in seq_len(size_step))
  {
    temp[i+1] <- temp[i] - theta*(temp[i] - mu)*dt + sigma*sqrt(dt)*rnorm(1)
  }
  
  sale <- rnorm(size_step+1, mean=20, sd=8)
  summary(lm(sale~temp))
  p_value[k] <- coef(summary(lm(sale~temp)))[2,4]
}
hist(p_value)


ほぼ一様な分布だ。さらに、 p値が5%を下回る割合を計算すると

> sum(p_value<0.05)/10^3
[1] 0.05

となったので、p値の割合(信頼度としての意味合いとしてのp値) としてもよさそうだ。

経路の可視化

ちなみに途中の1パスを可視化すると、以下のような感じ。

#plot
library("ggplot2")
data.frame(
  date=rep(Sys.Date()+seq_len(size_step+1), 2), 
  group=rep(c("temp", "sale"), each=size_step+1), 
  value=c(temp, sale)
) %>%
ggplot(aes(x=date, y=value, group=group, color=group)) + geom_line(size=2)

売上げ(平均回帰過程)を気温(平均回帰過程)で回帰してみる

業種によってはそうなのかもしれないが、売り上げ自体が平均回帰的な会社もあるかもしれない。そんなケースもあるはずだと思って同様の計算をしてみる。

size_step <- 100
sigma <- 5
dt <- 0.1
theta <- 1
mu_x <- 25
mu_y <- 35
p_value <- numeric(10^3)
for(k in 1:10^3)
{
  x <- rep(mu_x, size_step+1)
  y <- rep(mu_y, size_step+1)
  for(i in seq_len(size_step))
  {
    x[i+1] <- x[i] - theta*(x[i] - mu_x)*dt + sigma*sqrt(dt)*rnorm(1)
    y[i+1] <- y[i] - theta*(y[i] - mu_y)*dt + sigma*sqrt(dt)*rnorm(1)
  }
  
  p_value[k] <- coef(summary(lm(y~x)))[2,4]
}
hist(p_value, 100)

・・・明らかに有意となるケースが量産されている、これはやっちゃだめなやつだ!!!

> sum(p_value<0.05)/10^3
[1] 0.507
経路の可視化


見せかけの回帰の例

さらにちなみに、ダメなケース(ランダムウォーク同士の回帰、みせかけの回帰)だと、シードを固定しなくても、2〜3回も下記コードをまわせば即"みせかけの回帰"として有意な結果(xの回帰係数が有意)を作ることができる。

x <- cumsum(rnorm(100))
y <- cumsum(rnorm(100))

例えば以下で、xの係数が有意となっているが、これは明らかにみせかけの回帰の例だ。

> summary(lm(y~x))

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.3241 -2.2174  0.0802  2.5778  5.8766 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.53328    0.67048   6.761 9.94e-10 ***
x           -0.28741    0.05257  -5.467 3.48e-07 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.052 on 98 degrees of freedom
Multiple R-squared:  0.2337,	Adjusted R-squared:  0.2259 
F-statistic: 29.89 on 1 and 98 DF,  p-value: 3.484e-07

まとめ

解析的に計算している論文(Phillips(1986))だとランダムウォーク(正規分布を累積していくもの)同士の回帰がダメだと主張していたはずなので、それに該当しないケースはOK,あるいは個別に考えて対応となるのだろうか。

ここでは

  • 1: 正規分布 v.s. 平均回帰過程
  • 2: 平均回帰過程 v.s. 平均回帰過程

で回帰したが、2のケースは明らかにみせかけの回帰であると疑ってもよさそうだ。なので、過程(確率過程)に対するノイズが正規分布の累積(ランダムウォーク)になってる過程同士の回帰がまずそうだといえそうだ。

特に今回の売上のケースはよくよく考えると平均回帰過程でモデリングするのは筋が悪く、1のように単純に"ある想定される売り上げの水準からちょっとずれる"というような考え方をベースに持ってくるべきかなと*1

もっと簡単にいうと、気温の過程は平均回帰過程であるとした上では、売り上げをどうモデリングするか”が肝っぽい。

  • 1:売り上げ(今期) = (数値)×売り上げ(前期) + 正規分布 → アウト
  • 2:売り上げ(今期) = ある値(平均的な売り上げ) + 正規分布 → OK

ってかんじだ。1のケースが見せかけの回帰に相当するんで、イケイケなべんちゃー(売り上げが前期比で130%UP(上記,数値に相当)!)みたいな会社だと、1がよさそうですが、その会社さんだとみせかけの回帰出せると思ってて、「あのイケイケベンチャーの売上は・・・なんと気温に依存していた!」みたいな話になるんでねーかな、みせかけだけど。

*1:数式で考えると売上げに対するノイズ(バラツキ)が次の期の売上へと影響を与える形になるので筋が悪いなって

標準誤差と(ノンパラメトリック)ブートストラップのイケない関係

標準誤差と(ノンパラメトリック)ブートストラップについて、特に理論的なお話はおいておいて、直感的にどういうことなのかということをメモりたい。
それぞれ

  • (回帰分析の係数の)標準誤差: ある確率分布からのサンプルデータに対し推定されたある統計量(ここでは回帰係数)の、その(我々には見ることができない確率分布が持つ真の)統計量からのバラツキの大きさ
  • (ノンパラメトリック)ブートストラップ標準誤差: ブートストラップサンプルで推定されたある統計量(ここでは回帰係数)の、そのバラツキの大きさ

なので、これらは直感的に近しい値になるだろうということが期待されるわけです。これをやってみる。

まず、ふつーに回帰して標準誤差を取得すると・・・

> standard_error <- summary(lm(dist~speed, cars))$coef[2, 2]
> standard_error
[1] 0.4155128

となる。

次にノンパラメトリックブートストラップで標準誤差を計算してみる。コードは以下。

set.seed(71)
#ブートストラップのサイズを設定&結果格納用ベクトル定義
size_loop <- 10^4
coefs <- numeric(size_loop)
for(i in seq_len(size_loop))
{
  #サンプリングを実行
  bootstrap_index <- sample(nrow(cars), nrow(cars), replace=TRUE)
  #取得したサンプルに対してSpeed変数に対する回帰係数を取得
  coefs[i] <- summary(lm(dist~speed, cars[bootstrap_index, ]))$coef[2, 1]
}

その回帰係数の標準偏差がその係数のバラツキを表現するので計算してみると

> sd(coefs)
[1] 0.4086198

となり、直観どおり、大体、lm関数の返す標準誤差とあってるぞと。