R で米国 ETF パフォーマンス分析

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

library(tidyverse)

## 株価の取得
library(tidyquant)

## data.frame から xts 含む時系列データへの変換 (またはその逆)
library(timetk)

## パフォーマンス指標の計算
library(PerformanceAnalytics)

## ggplot2 関連
library(GGally)
library(ggrepel)
library(patchwork)

## データの表示の体裁を整えるパッケージ
library(formattable)

1

load("us-etf-research.RData")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

etfs <- tribble(
  ~class,        ~symbol,  ~index,
  "Equity",      "SPY",    "S&P 500",
  "Equity",      "IJH",    "S&P MidCap 400",
  "Equity",      "IJR",    "S&P SmallCap 600",
  "Equity",      "IWB",    "Russel 1000",
  "Equity",      "IWM",    "Russel 2000",
  "Equity",      "QQQ",    "NASDAQ 100",
  "Equity",      "DIA",    "Dow 30",
  "US Sector",   "XLE",    "Energy",
  "US Sector",   "XLF",    "Financials",
  "US Sector",   "XLU",    "Utilities",
  "US Sector",   "XLI",    "Industrial",
  "US Sector",   "XLK",    "Technology",
  "US Sector",   "XLV",    "Health",
  "US Sector",   "XLY",    "Consumer Discretionary",
  "US Sector",   "XLP",    "Consumer Staples",
  "US Sector",   "XLB",    "Materials",
  "Treasury",    "SHV",    "Short Term",
  "Treasury",    "SHY",    "1-3 Year",
  "Treasury",    "IEF",    "7-10 Year",
  "Treasury",    "TLT",    "20+ Year",
  "Treasury",    "TIP",    "Inflation Protected",
  "Corp Bond",   "LQD",    "Investment Grade",
  "Corp Bond",   "HYG",    "High Yield",
  "Corp Bond",   "PFF",    "Preffered Stock",
  "Real Estate", "IYR",    "Dow Jones U.S. Real Estate",
  "Commodity",   "GLD",    "Gold",
  "Commodity",   "SLV",    "Silver",
  "Commodity",   "USO",    "Crude Oil",
  "Forex",       "UUP",    "USD",
  "Forex",       "FXY",    "JPY",
  "Forex",       "FXE",    "EUR",
  "Forex",       "FXB",    "GBP",
  "Forex",       "FXF",    "CHF",
  "Volatility",  "VXX",    "S&P 500 VIX Short-Term Future"
)

1

rawdata <- tq_get(etfs$symbol, from = "1993-01-01")

1

head(rawdata) %>% mutate_if(is.numeric, round, digits = 2)

1
2
3
4
5

vxx <- read_csv("VXX.csv") %>% mutate(date = ymd(date))
rawdata <- rawdata %>%
  filter(!(symbol == "VXX" & date <= ymd("2019-01-29"))) %>%
  bind_rows(vxx) %>%
  arrange(symbol, date)

1
2
3
4
5
6

data <- rawdata %>%
  group_by(symbol) %>%
  tq_transmute(adjusted, mutate_fun = dailyReturn, col_rename = "ret") %>%
  slice(-1) %>% # リターン計算による先頭 0 を削除
  ungroup() %>%
  left_join(etfs, by = "symbol") # インデックスでもアクセスできるように結合しておく

1
2
3

data %>%
  group_by(symbol, index) %>%
  summarize(start = min(date), end = max(date), na = sum(is.na(ret)))

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

## 複数の銘柄の中で、最も遅い取引開始日を返す
get_start_date <- function(data) {
  data %>%
    group_by(symbol) %>%
    summarise(start_date =  min(date)) %>%
    pull(start_date) %>%
    max()
}

## get_start_date() の日付に全ての銘柄の開始日を合わせる
align_start_date <- function(data) {
  start_date <- get_start_date(data)
  filter(data, date >= start_date)
}

## 累積リターンとドローダウンを追加する
add_cumret_dd <- function(data) {
  ## PerformanceAnalytics::Drawdowns() 関数の warning を抑制
  suppressWarnings(
    data %>%
      group_by(symbol) %>%
      mutate(cum_ret = cumprod(1 + ret) - 1,
             drawdown = Drawdowns(ret)) %>%
      ungroup()
  )
}

## data.frame から xts へ変換
convert_to_xts <- function(data, name_col, value_col) {
  name_col <- enquo(name_col)
  value_col <- enquo(value_col)
  data %>%
    select(date, !!name_col, !!value_col) %>%
    pivot_wider(names_from = !!name_col, values_from = !!value_col) %>%
    tk_xts(-date, date)
}

## 凡例なしのシンプルな ggplot2 テーマ
theme_my <- function() {
  theme_light() + theme(
    plot.title = element_text(face = "bold.italic"),
    axis.title = element_blank(),
    legend.position = "none")
}

## 線グラフをプロット (累積リターン・ドローダウン・リターンで共通して利用)
plot_lines <- function(data, y_var, color_var, alpha = 0.6, offset = 0.2, label = TRUE) {
  y_var <- enquo(y_var)
  color_var <- enquo(color_var)

  min_date <- min(data$date)
  max_date <- max(data$date)
  range <- c(min_date, max_date + ((max_date - min_date) * offset))

  p <- data %>%
    group_by(!!color_var) %>%
    mutate(label = if_else(date == max(date), !!color_var, NA_character_)) %>%
    ggplot(aes(x = date, y = !!y_var, color = !!color_var)) +
    geom_line(alpha = alpha) +
    scale_x_date(limits = range) +
    theme_my()

  if (label) {
    p + geom_label_repel(aes(label = label), na.rm = TRUE, xlim = c(max_date, NA), size = 3)
  } else {
    p
  }
}

## {patchwork} を利用して複数の ggplot を 1 つに結合
plot_perf_summary <- function(data, alpha = 0.6, offset = 0.2) {
  p1 <- plot_lines(data, cum_ret, index, alpha, offset) +
    ggtitle("Cumulative Returns")

  p2 <- plot_lines(data, drawdown, index, alpha, offset, label = FALSE) +
    ggtitle("Drawdowns")

  p3 <- plot_lines(data, ret, index, alpha, offset, label = FALSE) +
    ggtitle("Daily Returns")

  p1 + p2 + p3 + plot_layout(ncol = 1, heights = c(2, 1, 1))
}

## {GGally} で散布図行列
plot_ggpairs <- function(data) {
  ## ベンチマークとの相関が高い順に並べる
  cors <- map_dfr(unique(data$symbol), function(symbol) {
    d1 <- filter(data, symbol == !!symbol) %>% pull(ret)
    d2 <- filter(data, symbol == !!benchmark) %>% pull(ret)
    cor <- cor(d1, d2)
    data.frame(symbol = symbol, cor = cor, stringsAsFactors = FALSE)
  })
  symbol_order <- arrange(cors, -cor) %>% pull(symbol)

  data %>%
    select(symbol, date, ret) %>%
    pivot_wider(names_from = symbol, values_from = ret) %>%
    select(-date) %>%
    ggpairs(aes(alpha = 0.1), columns = symbol_order, lower = list(continuous = "smooth")) +
    theme_my()
}

## 銘柄毎のパフォーマンス指標を計算して data.frame で返す
calc_perf_metrices <- function(data, etfs, benchmark = "SPY") {
  xts_data <- convert_to_xts(data, symbol, ret)
  rb <- xts_data[, benchmark]

  map_dfr(names(xts_data), function(symbol) {
    ra <- xts_data[, symbol]
    data.frame(
      Symbol  = symbol,
      Index   = etfs[etfs$symbol == symbol, ]$index,
      Sharpe  = round(as.numeric(SharpeRatio(ra, annualize = TRUE, FUN = "StdDev")), 2),
      Alpha   = round(CAPM.alpha(ra, rb), 5),
      Beta    = round(CAPM.beta(ra, rb), 2),
      Avg_Ret = as.numeric(Return.annualized(ra)),
      Cum_Ret = as.numeric(Return.cumulative(ra)),
      StdDev  = as.numeric(StdDev.annualized(ra)),
      MaxDD   = as.numeric(maxDrawdown(ra)),
      stringsAsFactors = FALSE
    )
  }) %>%
    mutate_at(vars(Avg_Ret, Cum_Ret, StdDev, MaxDD), percent, digits = 1) %>%
    arrange(-Sharpe)
}

1
2
3
4
5
6
7
8

asset_class <- list()
asset_class$symbols <- c("SPY", "IYR", "GLD", "TLT", "VXX")
asset_class$data <- data %>%
  filter(symbol %in% asset_class$symbols) %>%
  align_start_date()

asset_class$cor_plot <- plot_ggpairs(asset_class$data)
asset_class$cor_plot

1
2
3
4
5

equity <- list()
equity$data <- data %>%
  filter(class == "Equity") %>%
  align_start_date() %>%
  add_cumret_dd()

1
2

equity$perf_plot <- plot_perf_summary(equity$data, offset = 0.15)
equity$perf_plot

1
2

equity$perf_table <- calc_perf_metrices(equity$data, etfs, benchmark = "SPY")
equity$perf_table

1
2

equity$cor_plot <- plot_ggpairs(equity$data)
equity$cor_plot

1
2
3
4
5

us_sector <- list()
us_sector$data <- data %>%
  filter(class == "US Sector" | symbol == "SPY") %>% # ベンチマークとして SPY を含める
  align_start_date() %>%
  add_cumret_dd()

1
2

us_sector$perf_plot <- plot_perf_summary(us_sector$data, offset = 0.15)
us_sector$perf_plot

1
2

us_sector$perf_table <- calc_perf_metrices(us_sector$data, etfs, benchmark = "SPY")
us_sector$perf_table

1
2

us_sector$cor_plot <- plot_ggpairs(us_sector$data)
us_sector$cor_plot

1
2
3
4
5

treasury <- list()
treasury$data <- data %>%
  filter(class == "Treasury" | symbol == "SPY") %>%
  align_start_date() %>%
  add_cumret_dd()

1
2

treasury$perf_plot <- plot_perf_summary(treasury$data, offset = 0.12)
treasury$perf_plot

1
2

treasury$perf_table <- calc_perf_metrices(treasury$data, etfs, benchmark = "SPY")
treasury$perf_table

1
2

treasury$cor_plot <- plot_ggpairs(treasury$data)
treasury$cor_plot

1
2
3
4
5

corp_bond <- list()
corp_bond$data <- data %>%
  filter(class == "Corp Bond" | symbol == "SPY") %>%
  align_start_date() %>%
  add_cumret_dd()

1
2

corp_bond$perf_plot <- plot_perf_summary(corp_bond$data, offset = 0.15)
corp_bond$perf_plot

1
2

corp_bond$perf_table <- calc_perf_metrices(corp_bond$data, etfs, benchmark = "SPY")
corp_bond$perf_table

1
2

corp_bond$cor_plot <- plot_ggpairs(corp_bond$data)
corp_bond$cor_plot

1
2
3
4
5

real_estate <- list()
real_estate$data <- data %>%
  filter(class == "Real Estate" | symbol == "SPY") %>%
  align_start_date() %>%
  add_cumret_dd()

1
2

real_estate$perf_plot <- plot_perf_summary(real_estate$data, offset = 0.2)
real_estate$perf_plot

1
2

real_estate$perf_table <- calc_perf_metrices(real_estate$data, etfs, benchmark = "SPY")
real_estate$perf_table

1
2

real_estate$cor_plot <- plot_ggpairs(real_estate$data)
real_estate$cor_plot

1
2
3
4
5

commodity <- list()
commodity$data <- data %>%
  filter(class == "Commodity" | symbol == "SPY") %>%
  align_start_date() %>%
  add_cumret_dd()

1
2

commodity$perf_plot <- plot_perf_summary(commodity$data, offset = 0.1)
commodity$perf_plot

1
2

commodity$perf_table <- calc_perf_metrices(commodity$data, etfs, benchmark = "SPY")
commodity$perf_table

1
2

commodity$cor_plot <- plot_ggpairs(commodity$data)
commodity$cor_plot

1
2
3
4
5

volatility <- list()
volatility$data <- data %>%
  filter(class == "Volatility" | symbol == "SPY") %>%
  align_start_date() %>%
  add_cumret_dd()

1
2

volatility$perf_plot <- plot_perf_summary(volatility$data, offset = 0.2)
volatility$perf_plot

1
2

volatility$perf_table <- calc_perf_metrices(volatility$data, etfs, benchmark = "SPY")
volatility$perf_table

1
2

volatility$cor_plot <- plot_ggpairs(volatility$data)
volatility$cor_plot

1
2
3
4
5

forex <- list()
forex$data <- data %>%
  filter(class == "Forex" | symbol == "SPY") %>%
  align_start_date() %>%
  add_cumret_dd()

1
2

forex$perf_plot <- plot_perf_summary(forex$data, offset = 0.1)
forex$perf_plot

1
2

forex$perf_table <- calc_perf_metrices(forex$data, etfs, benchmark = "SPY")
forex$perf_table

1
2

forex$cor_plot <- plot_ggpairs(forex$data)
forex$cor_plot