分割投資 vs 一括投資

例えば、現在5000万円を持っているとして、 A: 毎年250万円ずつ投資して20年かけて投資額（元本）を5000万円にする B: 半分の2500万を一括投資して運用どちらのほうがリスクが少ないか？

リスクの測定として、20年間運用した後の総資産額についての期待標準偏差を比較する。投資ポートフォリオの特性として、1年あたりのリターンが、年率期待リターンmu, 年率期待標準偏差 sigmaの正規分布に従うi.i.d系列とする。ポートフォリオリターンをBootstrap法により繰り返し発生させて比較を行う。

運用後の総資産額は指数的な分布になることが予測されるため、対数変換した値を比較する。また、最終資産額についての95%信頼区間（2.5パーセンタイル値, 97.5パーセンタイル値）も合わせて求める。

対数変換した後の期待値、SDおよびその比率も示す。

#ライブラリ読み込み
library(ggplot2)
library(ggsci)

投資シミュレーションの関数定義

Simulation <- function(mu, sigma, year, initialPrice){
  # mu: 年率の期待リターン（算術平均リターン）
  # sigma: 年率の期待標準偏差
  # year: 投資する年数
  # initialPrice: 開始時の資産額
  
  # Parameters
  N <- year
  
  # Return Generation
  yearlyReturn <- rnorm(N, mu, sigma) # 年率リターンを正規分布（N(mu, sigma)）に従うi.i.d系列として発生
  
  # 分割して投資する(A)
  totalAssetA <- rep(0, N) # 総資産額の格納
  cash <- rep(0, N)        # キャッシュ格納用
  investPrice <- rep(0, N) # 投資評価額の格納用
  
  cash[1] <- initialPrice # はじめのcashを設定
  contributeAmountA <- initialPrice/N # 1回あたり(毎年)の投資額
  
  # N年分シミュレーション
  for(i in 1:N){
    # 年初に投資
    investPrice[i] <- investPrice[i]+contributeAmountA
    cash[i] <- cash[i] - contributeAmountA
    
    # リターン計算
    investPrice[i] <- investPrice[i]*(yearlyReturn[i]+1)
    if(investPrice[i]<=0){investPrice[i] <- 0} # 0以下（破産）になったら0を代入
    totalAssetA[i] <- investPrice[i] + cash[i]
    if(i == N){break}
    
    # 繰越
    investPrice[i+1] <- investPrice[i]
    cash[i+1] <- cash[i]
  }
  
  # 半分を一括投資して運用し続ける(B)
  contributeAmountB <- initialPrice/2 # はじめに投資する額
  investPrice <- rep(0, N) # 投資評価額の格納用
  investPrice[1] <- contributeAmountB
  for(i in 1:N){
    investPrice[i] <- max(0, investPrice[i] * (1+yearlyReturn[i]))
    if(i==N){break}
    investPrice[i+1] <- investPrice[i]
  }
    totalAssetB <- investPrice[N] + (initialPrice - contributeAmountB)
  
  return(c(totalAssetA[N], totalAssetB)) # A, Bの最終資産額を出力
}

Simulation(0.09, 0.2, 20, 5000) # mu=0.09, sigma=0.2としてシミュレーションした結果

## [1] 4614.643 6227.518

シミュレーションの反復試行

Bootstrap法による反復試行は、1組の(mu, sigma)あたり100,000回とする。

repeatN <- 100000 #1組のパラメータあたりの試行数
resultA <- rep(0, repeatN) # Aの最終資産額（対数変換後）格納用
resultB <- rep(0, repeatN) # Bの最終資産額（対数変換後）格納用

# 以下のmu, sigmaの総当りの組み合わせで検定を行う。
m <- c(0.03, 0.06, 0.09, 0.12, 0.15) # mu(vector)
s <- c(0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5) # sigma(vector)

# 結果の格納用
sd_ratio_df <- data.frame(matrix(rep(0, length(m) * length(s)), nrow=length(s)))
mean_ratio_df <- data.frame(matrix(rep(0, length(m) * length(s)), nrow=length(s)))
mean_A_df <- data.frame(matrix(rep(0, length(m) * length(s)), nrow=length(s)))
mean_B_df <- data.frame(matrix(rep(0, length(m) * length(s)), nrow=length(s)))
median_A_df <- data.frame(matrix(rep(0, length(m) * length(s)), nrow=length(s)))
median_B_df <- data.frame(matrix(rep(0, length(m) * length(s)), nrow=length(s)))
sd_A_df <- data.frame(matrix(rep(0, length(m) * length(s)), nrow=length(s)))
sd_B_df <- data.frame(matrix(rep(0, length(m) * length(s)), nrow=length(s)))
lCI_A_df <- data.frame(matrix(rep(0, length(m) * length(s)), nrow=length(s)))
uCI_A_df <- data.frame(matrix(rep(0, length(m) * length(s)), nrow=length(s)))
lCI_B_df <- data.frame(matrix(rep(0, length(m) * length(s)), nrow=length(s)))
uCI_B_df <- data.frame(matrix(rep(0, length(m) * length(s)), nrow=length(s)))
results <- data.frame(matrix(rep(NA, 4), nrow=1))[numeric(0), ]
colnames(results) <- c("method", "mu", "sigma", "value")

# ループ
for(k in 1:length(m)){ # muについてのループ
for(j in 1:length(s)){ # sigmaについてのループ
  for(i in 1:repeatN){ # bootstrap
    # 特定の(mu, sigma)についてのシミュレーション
    res <- Simulation(m[k], s[j], 20, 50000000) # 1試行での結果（最終資産額）
    resultA[i] <- log10(max(res[1],1)) # A: 積立投資
    resultB[i] <- log10(max(res[2],1)) # B: 一括投資
  }
  
  # 結果（対数変換した最終総資産額の統計量を求める）
  mean_A      <- mean(resultA)
  mean_B      <- mean(resultB)
  median_A    <- median(resultA)
  median_B    <- median(resultB)
  sd_A        <- sd(resultA)
  sd_B        <- sd(resultB)
  sd_ratio    <- (sd_B/mean_B) / (sd_A/mean_A)
  mean_ratio  <- mean_B - mean_A
  CI_A <- quantile(resultA, c(0.025, 0.975))
  CI_B <- quantile(resultB, c(0.025, 0.975))
  results <- rbind(results, 
                 data.frame(method=rep("Dollar Cost",repeatN), 
                            mu=rep(m[k], repeatN), 
                            sigma=rep(s[j], repeatN), 
                            value=resultA),
                 data.frame(method=rep("Bulk",repeatN), 
                            mu=rep(m[k], repeatN), 
                            sigma=rep(s[j], repeatN), 
                            value=resultB))

  # テーブルにまとめる
  mean_A_df[j,k]     <- mean_A
  mean_B_df[j,k]     <- mean_B
  median_A_df[j,k]   <- median_A
  median_B_df[j,k]   <- median_B
  sd_A_df[j,k]       <- sd_A
  sd_B_df[j,k]       <- sd_B
  sd_ratio_df[j,k]   <- sd_ratio
  mean_ratio_df[j,k] <- mean_ratio
  lCI_A_df[j,k]       <- CI_A[1]
  uCI_A_df[j,k]       <- CI_A[2]
  lCI_B_df[j,k]       <- CI_B[1]
  uCI_B_df[j,k]       <- CI_B[2]
}
}

library(ggplot2)
library(ggsci)
# 箱ひげ図でそれぞれのsigmaについての中央値、四分位を示す。
plot_est_ci_bp <- function(results, sigma){
  title = paste("Final wealth distribution after 20 years in sigma = ", sigma, sep="")
  df_plot <- results[results$sigma==sigma,]
  g <- ggplot(df_plot, aes(x=as.factor(mu), y=value, fill=method))+
    geom_boxplot() +
    scale_color_nejm() +
    xlab("mu") + ylab("Wealth") +
    labs(title=title) +
    ylim(c(6,NA))
  return(g)
}

plot_est_ci_bp(results, 0.05)

plot_est_ci_bp(results, 0.1)

plot_est_ci_bp(results, 0.2)

plot_est_ci_bp(results, 0.3)

## Warning: Removed 135 rows containing non-finite values (stat_boxplot).

plot_est_ci_bp(results, 0.4)

## Warning: Removed 2227 rows containing non-finite values (stat_boxplot).

plot_est_ci_bp(results, 0.5)

## Warning: Removed 9816 rows containing non-finite values (stat_boxplot).

結果のテーブルをプロット

library(ggplot2)
library(reshape2)
library(ggsci)

plot_result <- function(df, colname, rowname, ylabel, title_char){
  colnames(df) <- colname
  rownames(df) <- rowname
  df <- stack(df)
  df$sigma <- rep(s, length(m))
  df$mu <- df$ind
  g <- ggplot(df, aes(x=sigma, y=values, color=mu)) +
    geom_line() +
    scale_color_nejm() +
    ylab(ylabel) +
    labs(title=title_char)
  return(g)
}

plot_result2 <- function(df1, df2, colname, rowname, ylabel, title_char){
  colnames(df1) <- colname
  rownames(df1) <- rowname
  df1 <- stack(df1)
  df1$sigma <- rep(s, length(m))
  df1$mu <- df1$ind
  df1$Method <- "Dollar Cost"
  
  colnames(df2) <- colname
  rownames(df2) <- rowname
  df2 <- stack(df2)
  df2$sigma <- rep(s, length(m))
  df2$mu <- df2$ind
  df2$Method <- "Bulk"
  
  df <- rbind(df1, df2)
  
  g <- ggplot(df, aes(x=sigma, y=values, color=mu, linetype=Method)) +
    geom_line() +
    scale_color_nejm() +
    ylab(ylabel) +
    labs(title=title_char)
  return(g)
}

(g_sd <- plot_result(sd_ratio_df, m, s, "SD ratio: B/A", "SD ratio")) # meanで除した標準偏差の比率(B/A)

(g_mean <- plot_result(mean_ratio_df, m, s, "Mean ratio: B/A", "Mean ratio")) # 対数変換した総資産額のmeanの比率(B/A)

#(mean_A <- plot_result(mean_A_df, m, s, "Mean", "Mean of A (logged)")) # 対数化したA資産額のmean
#(mean_B <- plot_result(mean_B_df, m, s, "Mean", "Mean of B (logged)")) # 対数化したB資産額のmean
#(sd_A <- plot_result(sd_A_df, m, s, "SD", "StdDev of A (logged)")) # 対数化したA資産額の標準偏差
#(sd_B <- plot_result(sd_B_df, m, s, "SD", "StdDev of B (logged)")) # 対数化したB資産額の標準偏差
(plot_result2(mean_A_df, mean_B_df, m, s, "Wealth", "Mean Wealth (logged)"))

(plot_result2(sd_A_df, sd_B_df, m, s, "SD", "SD (logged)"))

ドルコスト平均法vs一括投資

分割投資 vs 一括投資

投資シミュレーションの関数定義

シミュレーションの反復試行

結果のテーブルをプロット