ドルコスト平均法のご利益をRで理解する

掲題の件、そういう時あると思います。

結論

まあ、ちょっと考えれば自明なんだが、以下です。

  • ドルコスト平均法は平均的なリターンを押し下げる(儲かる投資なら!)効果があるので嬉しくはない
  • ドルコスト平均法は最終的な儲けのバラツキ(標準偏差)を押し下げる効果があるので、これは不確実性を削減出来ているという意味で嬉しい

状況と結果

投資戦略①

①全期間(250日間)において毎日一定金額(1円)を投資した場合の最終的な儲けとそのバラツキ

> performance(s1)
[1] 258.46619  30.96698

投資戦略②

②初日に全額(250円)を投資した場合の最終的な儲けとそのバラツキ

> performance(s2)
[1] 266.92645  53.44526

それぞれのシミュレーションを複数回実行した際の密度描画(①=黒線、②=赤線)。 初日に全額ブチ込んだ戦略のほうが標準偏差がデカイ(当たり前)。 f:id:teramonagi:20211016151327p:plain

Code

library("purrr")
# Simulate Investment
simulate <- function(strategy, return){
  x <- 0
  for(i in seq_len(T)){
    x <- x + strategy[i]
    x <- (1 + return[i]) * x
  }
  x
}
performance <- function(s){
  c(mean(s), sd(s))
}

# Invenstment horizon (days)
horizon <- 250
# Mean Growth rate (/250 means Annualy -> Daily)
mean <- 0.07 / 250
# Volatility (/sqrt(250) means Annualy -> Daily)
vol <- 0.2 / sqrt(250)
size <- 10^3
# Generate return series
xs <- purrr::map(seq_len(size), ~ rnorm(horizon, mean=mean, sd = vol))

# Check the performance of different strategies
s1 <- purrr::map_dbl(xs, ~ simulate(rep(1, horizon), .x))
s2 <- purrr::map_dbl(xs, ~ simulate(c(horizon, rep(0, horizon-1)), .x))
performance(s1)
performance(s2)

plot(density(s1), col=1)
lines(density(s2), col=2)

最も近い値を持つ要素のインデックスを返す

掲題の件、あると思います。

例えばこういうやつ(z(=3)の値に最も近いxの要素(=2.9)を探す)。

> z <- 3
> x <- c(1,2,2.9,3.5,4)
> index <- which.min(abs(z - x))
> x[index]
[1] 2.9

これをzがベクトルのときでも動くようにしたいと考えた2実装が以下。

1つ目は単純に sapply を使っていて、2つ目は findInterval区間を見つけつつちょいと補正(単純にやると最も近い値にならないときがある)をしたもの。

1億個要素を持つベクトル(x)と、複数個値を持ってる適当なベクトル(y)を使って2つのやり方の時間を測ると

> x <- sort(runif(10^8, min=1, max=4))
> ys <- seq(2,3.9,by=0.3)
> system.time({ ce_1 <- sapply(ys, function(z)which.min(abs(z - x))) })
   ユーザ   システム       経過  
     2.592      1.650      4.362 
> system.time({ ce_2 <- findInterval(ys, 0.5*(x[-length(x)] + x[-1])) + 1 })
   ユーザ   システム       経過  
     1.285      0.544      1.830 

findInterval で書いたほうが早いんで良さそう。

念の為結果をチェックすると

> all(ce_1 == ce_2)
[1] TRUE

ちゃんとあってる。

まぁ、パッケージにするまでもないし、自分で関数として持っておくか〜という決断をする。

大きめのSparse行列をnormalizeしたい時の書き方

掲題の件、あると思います。

私は今まで行列っぽいデータを normalize するとき*1は以下のように書いていました。

最後の結果を見たらわかるように行ごとに和をとると1になるようになっています。

> x <- matrix(1:4,2)
> x
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> x/rowSums(x)
          [,1]      [,2]
[1,] 0.2500000 0.7500000
[2,] 0.3333333 0.6666667
> x <- matrix(1:4,2)
> x
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> x/rowSums(x)
          [,1]      [,2]
[1,] 0.2500000 0.7500000
[2,] 0.3333333 0.6666667

適当に作ったSparse行列でもこれで動くと言えば動いている。

> library("Matrix")
> x <- as(matrix(c(1,2,0,4,0,0,7,8,9),3),"sparseMatrix")
> x
3 x 3 sparse Matrix of class "dgCMatrix"
          
[1,] 1 4 7
[2,] 2 . 8
[3,] . . 9
> x/Matrix::rowSums(x)
3 x 3 sparse Matrix of class "dgCMatrix"
                                   
[1,] 0.08333333 0.3333333 0.5833333
[2,] 0.20000000 .         0.8000000
[3,] .          .         1.0000000

でも、行列のサイズがでかくなるとアウト、なんか一回Denseな行列にしてるんだろうと思われる。

ここの例は恣意的 データに見えるがもともと俺が使ってハマったデータがこんなんだったのでこうなっている。

> size <- 10^5
> i <- c(sample(1:size, 3*10^7, replace=TRUE), size)
> j <- c(sample(1:size, 3*10^7, replace=TRUE), size)
> x <- sparseMatrix(i,j,x=rnorm(3*10^7 + 1), dims=rep(size + 1, 2))
> x <- as(x + t(x), "symmetricMatrix")
> str(x)
Formal class 'dsCMatrix' [package "Matrix"] with 7 slots
  ..@ i       : int [1:29910604] 2 23 14 6 13 28 18 24 0 19 ...
  ..@ p       : int [1:100002] 0 0 0 0 0 0 0 0 0 0 ...
  ..@ Dim     : int [1:2] 100001 100001
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num [1:29910604] -1.511 -0.878 0.777 1.258 1.185 ...
  ..@ uplo    : chr "U"
  ..@ factors : list()
> x/Matrix::rowSums(x)
 エラー:  ベクトルのメモリを使い切りました (上限に達した?) 

こういうときは等価な計算のこれを使おう、これならいける

> Matrix::Diagonal(x=1 / Matrix::rowSums(x)) %*% x
100001 x 100001 sparse Matrix of class "dgCMatrix"
                                                                                                                              
 [1,] . . . . . . . . . . . . . . . . . . . . . .  .          . . . . . . . . . . . . . . . . .  .        . . . . . . . ......
 [2,] . . . . . . . . . . . . . . . . . . . . . .  .          . . . . . . . . . . . . . . . . .  .        . . . . . . . ......
 [3,] . . . . . . . . . . . . . . . . . . . . . . -0.08169018 . . . . . . . . . . . . . . . . .  .        . . . . . . . ......

ちゃんと計算も合う(確認)

> x <- as(matrix(c(1,2,0,4,0,0,7,8,9),3),"sparseMatrix")
> x/Matrix::rowSums(x)
3 x 3 sparse Matrix of class "dgCMatrix"
                                   
[1,] 0.08333333 0.3333333 0.5833333
[2,] 0.20000000 .         0.8000000
[3,] .          .         1.0000000
> Matrix::Diagonal(x=1 / Matrix::rowSums(x)) %*% x
3 x 3 sparse Matrix of class "dgCMatrix"
                                   
[1,] 0.08333333 0.3333333 0.5833333
[2,] 0.20000000 .         0.8000000
[3,] .          .         1.0000000

元ネタ

*1:各行ごとに足したら1になる、でnormalizeを定義。他L2ノルムなどによるnormalizeも同様

Pythonの純粋仮想関数は引数の数適当に変えて実装してもOK

掲題の件、そういいうことです。 例えば適当に get() メソッドの引数を追加して実装しても

from abc import ABCMeta, abstractmethod

class Hoge(metaclass=ABCMeta):
    @abstractmethod
    def get(self):
        pass

class Moge(Hoge):
    def get(self, x = None, y = None):
        print("hi")

ちゃんと動く&元のクラスのインスタンスは作れない

>>> x = Moge()
>>> x.get(y=123)
hi
>>> 
>>> Hoge()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: Can't instantiate abstract class Hoge with abstract methods get

米国株式とUSDJPYの2020年1月1日からの相関構造

f:id:teramonagi:20200820065748p:plain
米国株式とUSDJPYの2020年1月1日からの相関構造

ドル円(USDJPY)が 0.6% で ダウ平均が2.8% くらいのリスク(標準偏差、日率)

> 100*sd(df$`usdjpy=x`)
[1] 0.6186612
> 100*sd(df$dji)
[1] 2.842537

全体

library("dplyr")
library("tidyr")
library("stringr")
library("tidyquant")
library("PerformanceAnalytics")
get_from_yahoo <- function(x, range_date){
  # Why max(range_date) +1? because: from <= x < to
  df <- tq_get(x, from = min(range_date), to = max(range_date) + 1)
  dplyr::select(df, symbol, date, adjusted) %>%
    rename(price=adjusted) %>%
    mutate(symbol = str_to_lower(str_replace_all(symbol, "\\^", ""))) %>%
    arrange(date)
}
arithmetic_return <- function(x){
  c(NA, x[-1]/x[-length(x)] - 1)
}
# https://finance.yahoo.com/quote/%5EDJI/
range_date <- c(lubridate::ymd("20200101"), lubridate::ymd("20200818"))
symbols <- c("^DJI", "USDJPY=X")
df <- purrr::map_df(symbols, ~ get_from_yahoo(.x, range_date)) %>%
  tidyr::pivot_wider(names_from = symbol, values_from = price) %>%
  arrange(date) %>%
  mutate_at(
    vars(-"date"),
    arithmetic_return
  ) %>%
  na.omit()

chart.Correlation(dplyr::select(df, -date), histogram=TRUE, pch=19)
100*sd(df$`usdjpy=x`)
100*sd(df$dji)

米国株式市場の2020年1月1日からの相関構造

f:id:teramonagi:20200819124957p:plain
米国株式市場の2020年1月1日からの相関構造

すげぇ強かった(おしまい)

  • "^DJI" : ダウ平均
  • "^RUT": ラッセル指数
  • "^GSPC": S&P500
  • "^IXIC": ナスダック指数

再現用のCode

library("dplyr")
library("tidyr")
library("stringr")
library("tidyquant")
library("PerformanceAnalytics")
get_from_yahoo <- function(x, range_date){
  # Why max(range_date) +1? because: from <= x < to
  df <- tq_get(x, from = min(range_date), to = max(range_date) + 1)
  dplyr::select(df, symbol, date, adjusted) %>%
    rename(price=adjusted) %>%
    mutate(symbol = str_to_lower(str_replace_all(symbol, "\\^", ""))) %>%
    arrange(date)
}
arithmetic_return <- function(x){
  c(NA, x[-1]/x[-length(x)] - 1)
}
# https://finance.yahoo.com/quote/%5EDJI/
range_date <- c(lubridate::ymd("20200101"), lubridate::ymd("20200818"))
symbols <- c("^DJI", "^RUT", "^GSPC", "^IXIC")
df <- purrr::map_df(symbols, ~ get_from_yahoo(.x, range_date)) %>%
  tidyr::pivot_wider(names_from = symbol, values_from = price) %>%
  arrange(date) %>%
  mutate_at(
    vars(-"date"),
    arithmetic_return
  ) %>%
  na.omit()

chart.Correlation(dplyr::select(df, -date), histogram=TRUE, pch=19)

特定の列だけ外して処理したい場合には . (dot) を名前に持つ列と lsを組み合わせると良い

株式会社ホクソエムの社長から教えていただいた。 いちいち dplyr::select なぞせんでもこうするだけで .y を除いた処理を実行できる。

> df <- data.frame(
+   x = 1:3,
+   .y = 2:4
+ )
> df
  x .y
1 1  2
2 2  3
3 3  4
> ls(df)
[1] "x"