How to backtest? : Sector rotation example

[ 참고자료: 섹터 ETF를 이용한 섹터 로테이션 전략 ]

llibrary(PerformanceAnalytics)
library(quantmod)

symbols = c("091230.KS", "091220.KS", "098560.KS", "139220.KS", "139290.KS", "139270.KS", "139250.KS", "139230.KS", "139240.KS", "139280.KS", "139260.KS", "143860.KS", "152180.KS", "157490.KS",
"157510.KS", "157500.KS", "157520.KS", "227550.KS", "227560.KS", "227540.KS", "140710.KS", "140700.KS", "148070.KS")

getSymbols(symbols, src='yahoo', from = '2011-01-01', to = Sys.Date() )

먼저, PerformanceAnalytics 와 quantmod 패키지를 불러옵니다.
만약 설치 되지 않은 패키지라면
install.packages("패키지이름") 으로 설치 하시면 됩니다.

다음은 symbols 공간에 다운 받을 주식 종목들의 티커를 입력합니다. R 패키지에서는 크게 야후와 구글 파이낸스를 통해 주가 데이터를 받을 수 있습니다. (개인적으로 수정주가를 받을 수 있는 야후를 더 선호합니다.) 국내 거래소 기준 TICKER 6자리 & .KS 를 추가해주면 됩니다. 다른 주식으로 바꾸시고 싶으면 티커만 확인한 후 저 포맷으로 바꿔주시면 됩니다.

기존 연구에 사용했던 ETF 들을 그대로 입력하였습니다.

그 후, getSymbols 명령어를 통해 야후에서 데이터를 받아옵니다. Sys.Date() 는 컴퓨터 기준 시간, 즉 오늘까지 데이터를 받는 다는 뜻입니다.

prices = list()
for(i in 1:length(symbols)) {
  tmp = Ad(get(symbols[[i]]))
  prices[[i]] = tmp
}
prices = do.call(cbind, prices)

다음 받은 주가들을 prices 라는 매트릭스에 집어 넣는 과정입니다.
Ad, 즉 수정주가를 받아 옵니다.

ret_d = as.matrix(Return.calculate(prices))
ret_m = apply.monthly(ret_d, Return.cumulative)
ret_m[which(ret_m==0)] = NA

ret_sec = ret_m[, 1: (ncol(ret_m) -1 )]
ret_cash = ret_m[, ncol(ret_m)]

먼저 Return.calculate 명령어는 수익률을 계산하는 명령어 입니다.
기존에 사용했던 Delt는 하나의 column 밖에 계산되지 않는 반면, 이 명령어는 matrix 형태의 데이터를 산술 혹은 로그 수익률로 계산이 가능합니다.

20년, 2,000 종목의 daily return 을 계산하는데도 몇초면 될 정도로, 속도 또한 굉장합니다.
단, 해당 명령어를 쓰면 xts 형태로 저장이 되어, 나중에 계산을 하는데 번거로워 집니다. as.matrix는 괄호 안의 내용을 매트릭스 형태로 바꿔준다는 명령어 입니다.

apply.monthly는 괄호 안의 매트릭스를 월간 수익률로 바꿔주는 기능입니다. apply.yearly 를 사용하시면 년간 수익률로 자동으로 바꿔줍니다.

which(ret_m) == 0 은, ret_m 중에 0인 데이터를 찾는 명령어입니다. 아직 상장되지 않은 종목의 경우 daily return은 없는 반면, apply.monthly 명령어를 통해 월간 수익률 형태로 바꿔주면 0 으로 표시 되어 버립니다. 없는 종목인데 있는 종목 처럼 보이는 것입니다. 그래서 다시 이 수익률을 NA (Not available) 로 바꿔줍니다.

ret_sec에는 ret_m (월간수익률) 의 첫째부터 (마지막-1) 까지 열의 데이터를 가지고 옵니다. 마지막 열 데이터는 Cash 대용인 국고채(TICKER: 140700) 수익률 이기 때문입니다.

period = 12

ret_sec_mom =  matrix(0,nrow(ret_sec),ncol(ret_sec))
  for ( i in period : nrow(ret_sec) ) {
    ret_sec_mom[i, ] = apply( (ret_sec [(i-(period-1)):i , ] + 1), 2, prod) - 1
  }

먼저 셋팅하고 싶은 모멘텀 기간을 정합니다. (여기서는 12개월로 정하였습니다.)
matrix 함수를 통해 zero 매트릭스를 먼저 설정해 줍니다.

ret_sec [(i - (period-1)) : i, ] 는 ret_sec 중 현재 부터 (period-1) 까지의 행, 즉 현재가 15년 12월 이라면 15년 12월부터 11행 전인 15년 1월 까지의 수익률을 가져옵니다. 그 데이터에 각각 1 을 더해줍니다. 그 후 2 를 해주는 이유는 apply 명령어에서 1은 행 단위,  2는 열 단위 계산을 행합니다. 우리는 Time Series 형태로 데이터를 가지고 있기 때문에 열단위로 계산을 해 주어야 하므로 2를 입력합니다. prod는 이 값들을 곱하라는 것입니다. 이 곱한 값에  -1 을 빼줍니다.

간단하게 생각하면 이렇습니다. 3개월 간 수익률이 각각 +1%, +2%, +3% 라고 생각해 보죠. 누적 수익률을 구하는 방법은 { (1+0.01) * (1+0.02) * (1+0.03) } - 1 = 0.0611 입니다. 이것을 프로그래밍으로 나타낸 것이 위의 형태입니다. 엑셀로 먼저 해보시면 쉽게 이해가 되실 겁니다.
(더 간단하게 하는 법도 있겠지만, 여기까지가 제 머리의 한계이고 어째튼 돌아는 가니까....)

ret_port = matrix(0, nrow(ret_sec), 1)
rownames(ret_port) = rownames(ret_m)

for (i in period : (nrow(ret_m) - 1) ) {
  
  K = which( (ret_sec_mom[i,] > 0) & (rank(- ret_sec_mom[i, ]) <= 5) )
  
  eq = mean(ret_sec[i+1, K]) * length(K) / 5 
  eq[is.nan(eq)] = 0
  bd = (1 - length(K)/5) * ret_cash[i+1 ]
  
  ret_port[i+1] = eq + bd
  
}

ret_port = as.matrix( ret_port[ which(rownames(ret_port) == "2012-01-31") : nrow(ret_port), ] )
colnames(ret_port) = c("Sector Momentum")

먼저 계산된 값들이 저장될 ret_port 매트릭스를 만들어 줍니다. rownames 명령어를 통해 각 행이름애 date 를 저장해 줍니다.

ret_sec_mom[i,] > 0 은 누적 수익률이 0 보다 큰 종목 (Time series momentum),
rank ( -ret_sec_mom[i, ] ) <= 5 은 누적 수익률의 순위가 5보다 작은 종목 (Cross sectional momentum) 을 찾습니다. 말씀 드린듯이 R 에서 rank는 기본적으로 오름차순이고, 우리는 수익률이 높을 수록 좋기 때문에 - 값을 입력하여 내림차순의 형태로 계산합니다.

& 구문 (and) 을 통해 두 조건 모두 만족하는 종목들의 열 위치를 K 에 저장합니다.

length(K) 는 K 에 저장된 숫자의 총 갯수입니다. length(K) / 5 는 총 주식에 투자될 비중 이겠죠. K 에 저장된 종목이 5개 라면, 주식에 100% 가 투자될 것이며, 종목이 4 개 라면 80%가 투자될 것입니다. 나머지 20%는 cash 혹은 채권에 투자되겠죠.

ret_sec[i+1, K] 를 통해 선정된 종목의 다음월 수익률을 구합니다. mean 값을 취해주면 동일가중 포트폴리오의 수익률이 계산 되고, length(K) / 5 값까지 곱해주면, 주식 투자 비중 까지 고려된 동일가중 포트폴리오의 수익률이 계산됩니다.

K에 저장된 값이 없을 경우, length(K) / 5 가 NaN 값이 나오고,  eq에 NaN 값이 입력되어 계산이 불가하게 됩니다. 이런 경우를 위해 is.nan 함수를 통해 NaN 형태일 경우 eq 를 0 으로 만들어 줍니다. (주식 비중이 0 이었으니, 수익률도 0 이겠죠)

마지막으로 bd는 다음월 채권 투자 비중이 고려된 채권 수익률을 구합니다. 만약 주식에 전부 투자 했다면 (K=5), 1 - length(K)/5 는 0, 즉 채권 투자 비중은 0%가 될 것입니다.

마지막으로 ret_port[i+1, ] 에 주식과 채권 투자 수익률의 합을 저장해 줍니다.

charts.PerformanceSummary(ret_port)
apply.yearly(ret_port, Return.cumulative)
barplot(t(apply.yearly(ret_port, Return.cumulative)))

rbind(Return.cumulative(ret_port), table.AnnualizedReturns(ret_port), maxDrawdown(ret_port))

이제부터는 결과물을 확인하는 명령어 들입니다.
먼저, charts.PerformanceSummary는 누적수익률, 기간별 수익률, Drawdown을 함께 나타내 주는 그래프 입니다. rowname 이 date 형태가 아닐 경우, 오류가 납니다.

apply.yearly는 말씀 드린것 처럼 연간 수익률로 바꿔주는 명령어 입니다.

위에서 구한 값을 barplot을 통해 막대그래프로 표현합니다. t() 를 통해 연간 수익률을 transpose 해주어야 연도별 수익률이 나옵니다.

Return.cumulative는 누적 수익률을, table_AnnualizedReturns는 연환산 수익률, 변동성, 샤프지수를, maxDrawdown은 MDD를 구해줍니다. rbind함수를 통해 해당 값들을 행으로 묶어 줍니다.

colnames(ret_sec)[which( (ret_sec_mom[nrow(ret_sec_mom), ] > 0) & (rank(- ret_sec_mom[nrow(ret_sec_mom), ]) <= 5) ) ]

마지막으로, 현재 시점에서 투자할 종목들을 선택합니다. which 를 조건에 충족하는 종목들의 위치를 찾고, colnames 를 통해 조건에 만족하는 종목들의 TICKER 를 반환합니다. 해당 종목들에 20% 씩 투자하면 됩니다. (만일 5 종목이 안될 경우 나머지 비중만큼 채권에 투자하면 됩니다.)

위에 적힌 코드들을 쭉 쓰면 아래와 같습니다.
백테스트용 데이터를 인터넷에서 긁어 오므로,
패키지만 설치되어 있다면, 동일한 결과가 나올 것입니다.

library(PerformanceAnalytics)
library(quantmod)

symbols = c("091230.KS", "091220.KS", "098560.KS", "139220.KS", "139290.KS", "139270.KS", "139250.KS", "139230.KS", "139240.KS", "139280.KS", "139260.KS", "143860.KS", "152180.KS", "157490.KS",
"157510.KS", "157500.KS", "157520.KS", "227550.KS", "227560.KS", "227540.KS", "140710.KS", "140700.KS", "148070.KS")

getSymbols(symbols, src='yahoo', from = '2011-01-01', to = Sys.Date() )

prices = list()
for(i in 1:length(symbols)) {
  tmp = Ad(get(symbols[[i]]))
  prices[[i]] = tmp
}
prices = do.call(cbind, prices)

ret_d = as.matrix(Return.calculate(prices))
ret_m = apply.monthly(ret_d, Return.cumulative)
ret_m[which(ret_m==0)] = NA

ret_sec = ret_m[, 1: (ncol(ret_m) -1 )]
ret_cash = ret_m[, ncol(ret_m)]

period = 12

ret_sec_mom =  matrix(0,nrow(ret_sec),ncol(ret_sec))
  for ( i in period : nrow(ret_sec) ) {
    ret_sec_mom[i, ] = apply( (ret_sec [(i-(period-1)):i , ] + 1), 2, prod) - 1
  }


###############################

ret_port = matrix(0, nrow(ret_sec), 1)
rownames(ret_port) = rownames(ret_m)

for (i in period : (nrow(ret_m) - 1) ) {
  
  K = which( (ret_sec_mom[i,] > 0) & (rank(- ret_sec_mom[i, ]) <= 5) )
  
  eq = mean(ret_sec[i+1, K]) * length(K) / 5 
  eq[is.nan(eq)] = 0
  bd = (1 - length(K)/5) * ret_cash[i+1 ]
  
  ret_port[i+1] = eq + bd
  
}

ret_port = as.matrix( ret_port[ which(rownames(ret_port) == "2012-01-31") : nrow(ret_port), ] )
colnames(ret_port) = c("Sector Momentum")

###############################

charts.PerformanceSummary(ret_port)
apply.yearly(ret_port, Return.cumulative)
barplot(t(apply.yearly(ret_port, Return.cumulative)))

rbind(Return.cumulative(ret_port), table.AnnualizedReturns(ret_port), maxDrawdown(ret_port))

###############################

colnames(ret_sec)[which( (ret_sec_mom[nrow(ret_sec_mom), ] > 0) & (rank(- ret_sec_mom[nrow(ret_sec_mom), ]) <= 5) ) ]