library(quantmod) library(PerformanceAnalytics)
tickers <- c(“AAPL”, “MSFT”, “SBUX”) n <- length(tickers)
data <- new.env() getSymbols(tickers, src = ‘yahoo’, from = “2007-01-01”, to = “2018-12-31”, env = data, auto.assign = T) names(data) head(data$AAPL)
stock.3 <- c(data\(AAPL\)AAPL.Adjusted, data\(MSFT\)MSFT.Adjusted, data\(SBUX\)SBUX.Adjusted)
con = gzcon(url(‘http://www.systematicportfolio.com/sit.gz’,‘rb’)) source(con) close(con) # bt.prep(data,align = ‘remove.na’,fill.gaps = T) names(data) head(data\(prices) head(data\)weight)
prices <- data$prices ret <- prices / mlag(prices) - 1 head(ret) ret <- ret[-1, ]
ret.sample <- ret[‘2007/2009’, ] avg.ret <- colMeans(coredata(ret.sample)) # cov.d <- cov(coredata(ret.sample)) ones.vec <- matrix(rep(1,n), ncol = 1) # min var portfolio numerator <- solve(cov.d) %% ones.vec denominator <- t(ones.vec)%%numerator mvp.w <- numerator / as.numeric(denominator) mvp.w # use historical returns to compute average return mvp.ret <- sum(as.numeric(mvp.w)avg.ret) #compute mvp standard deviation library(tidyverse) mvp.std <- t(mvp.w)%%cov.d%*%mvp.w %>% as.numeric %>% sqrt
mvp.std
ret.m.dplyr <- prices %>% to.monthly(indexAt = “lastof”, OHLC = FALSE) %>% data.frame(date = index(.)) %>%
remove_rownames() %>%
gather(stock, price, -date) %>%
group_by(stock) %>%
mutate(returns = log(price) - log((lag(price)))) %>% select(-price) %>%
spread(stock, returns) %>%
slice(-1) %>% head(ret.m.dplyr)
ret.w.dplyr <- prices %>% to.weekly(indexAt = “lastof”, OHLC = FALSE) %>% data.frame(date = index(.)) %>% remove_rownames() %>%
gather(stock, price, -date) %>%
group_by(stock) %>% mutate(returns = log(price) - log((lag(price)))) %>% select(-price) %>% spread(stock, returns) %>% slice(-1) %>% head(ret.w.dplyr)
#min var monthly ret.sample1 <- ret.m.dplyr
avg.ret1 <- colMeans(coredata(ret.sample1)) cov.m <- cov(coredata(ret.sample1)) ones.vec <- matrix(rep(1,n), ncol = 1) #min var portfolio numerator1 <- solve(cov.m) %% ones.vec denominator1 <- t(ones.vec)%%numerator1 mvp.m1 <- numerator1 / as.numeric(denominator1) mvp.m1 # use historical returns to compute average return mvp.ret1 <- sum(as.numeric(mvp.m1)avg.ret1) #compute mvp standard deviation library(tidyverse) mvp.std1 <- t(mvp.m1)%%cov.m%*%mvp.m1 %>% as.numeric %>% sqrt
mvp.std1
#min var weekly ret.sample2 <- ret.w.dplyr avg.ret2 <- colMeans(coredata(ret.sample2)) # cov.w <- cov(coredata(ret.sample2)) ones.vec <- matrix(rep(1,n), ncol = 1) # min var portfolio numerator2 <- solve(cov.w) %% ones.vec denominator2 <- t(ones.vec)%%numerator2 mvp.wk <- numerator2 / as.numeric(denominator2) mvp.wk # use historical returns to compute average return mvp.ret2 <- sum(as.numeric(mvp.wk)avg.ret2) #compute mvp standard deviation library(tidyverse) mvp.std2 <- t(mvp.wk)%%cov.w%*%mvp.wk %>% as.numeric %>% sqrt
mvp.std2
#================================= mu<-0.06/12 return <- ret.m.dplyr Ax <- rbind(2cov(return), colMeans(return), rep(1, ncol(return))) Ax <- cbind(Ax, rbind(t(tail(Ax, 2)), matrix(0, 2, 2))) b0 <- c(rep(0, ncol(return)), mu, 1) out<-solve(Ax, b0) wgt<-out[1:3] wgt sum(wgt) ret.out<-sum(wgtcolMeans(return)) ret.out.annual<-ret.out12 ret.out.annual std.out<-sqrt(t(wgt)%%cov(return)%%wgt) std.out.annual<-std.outsqrt(12) std.out.annual #==================================================================================================== # Or write your own function to find min var for a specified return mu; # Reference: Introduction to R for Quantitative Finance: Chapter 2, p31 # Solve a diverse range of problems with R, one of the most powerful tools for quantitative finance #==================================================================================================== return = ret.m.dplyr #specified portfolio return: mu mu=0.06/12
minvariance <- function(return, mu) { #return <- log(tail(assets, -1) / head(assets, -1)) Ax <- rbind(2cov(return), colMeans(return), rep(1, ncol(return))) Ax <- cbind(Ax, rbind(t(tail(Ax, 2)), matrix(0, 2, 2))) b0 <- c(rep(0, ncol(return)), mu, 1) zx<-solve(Ax, b0) weight<-zx[1:ncol(return)] ret.out<-sum(weightcolMeans(return)) std.out<-sqrt(t(wgt)%%cov(return)%%wgt) list(weight=weight, rtn=ret.out, sd=std.out) }
minvariance(return, mu)