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)
#Homework2 #Q1 riskfreerate assumed 0.01
#no short sales constraint
Amat1 = cbind(rep(1,3),mu, diag(1,nrow=3)) # set the constraints matrix t(Amat1)
muP1 = seq(min(mu)+.0001,max(mu)-.0001,length=300)
library(quadprog) sdP1 = muP1 # set up storage for standard deviations of portfolio returns weights1 = matrix(0,nrow=300,ncol=3) # storage for portfolio weights
i=1 for (i in 1:length(muP1)) # find the optimal portfolios for each target expected return { bvec1 = c(1,muP1[i], rep(0,3)) # constraint vector result = solve.QP(Dmat=2*Sigma,dvec=rep(0,3),Amat=Amat1,bvec=bvec1,meq=0) sdP1[i] = sqrt(result\(value) weights1[i,] = result\)solution }
pdf(return,width=6,height=5)
plot(sdP1,muP1,type=“l”,xlim=c(0,0.25),ylim=c(0,0.08),lty=1) # plot
#the efficient frontier (and inefficient frontier) par(new=TRUE) plot(sdP,muP,type=“l”,xlim=c(0,0.25),ylim=c(0,0.08),lty=1, col=“green”) # plot mufree = 0.005 # input value of risk-free interest rate points(0,mufree,cex=3,pch="") # show risk-free asset sharpe =( muP-mufree)/sdP # compute Sharpe ratios ind = (sharpe == max(sharpe)) # Find maximum Sharpe ratio options(digits=3) weights[ind,] # Find tangency portfolio# show line of optimal portfolios lines(c(0,2),mufree+c(0,2)(muP[ind]-mufree)/sdP[ind],lwd=2,lty=3)
con = gzcon(url(‘https://github.com/systematicinvestor/SIT/raw/master/sit.gz’, ‘rb’)) source(con) close(con)
n <- dim(return)[2] constraints = new.constraints(n, lb = -Inf, ub = +Inf)
constraints = add.constraints(rep(1, n), 1, type = ‘=’, constraints)
ia <- create.historical.ia(return, 365)
weight <- min.risk.portfolio(ia, constraints)
ia <- create.historical.ia(return, 365)
constraints = new.constraints(n, lb = 0, ub = 1) constraints = add.constraints(diag(n), type=‘>=’, b=0, constraints) constraints = add.constraints(diag(n), type=‘<=’, b=1, constraints)
constraints = add.constraints(rep(1, n), 1, type = ‘=’, constraints) # weight <- min.risk.portfolio(ia, constraints)
library(corpcor) library(lpSolve) ifelse(!require(corpcor), install.packages(“corpcor”), library(corpcor)) ifelse(!require(lpSolve), install.packages(“lpSolve”), library(lpSolve)) ef = portopt(ia, constraints, 1, ‘Efficient Frontier’) ef
library(pacman) p_load(timeSeries, PerformanceAnalytics, fPortfolio) return = return # monthly ret.ts_m <- timeSeries(ret.m.dplyr,) chart.CumReturns(ret.ts_m, legend.loc = ‘topleft’, main = ’’)
plot(portfolioFrontier(ret.ts_m))
Spec = portfolioSpec() setSolver(Spec) = “solveRshortExact” setTargetReturn(Spec) = mean(colMeans(ret.ts_m))
Data_m<-ret.ts_m nAssets_m = getNAssets(portfolioData(Data_m)) Weights_m <- rep(1/nAssets_m, times = nAssets_m) covRisk(Data_m, Weights_m) varRisk(Data_m, Weights_m, alpha = 0.05)###VaR of equal weights portfolio cvarRisk(Data_m, Weights_m, alpha = 0.05)### CVaR for equal weight portfolio
return_m = ret.m.dplyr setRiskFreeRate(Spec) <- 0.01/12 tgPortfolio <- tangencyPortfolio(ret.ts, Spec, constraints = “Short”) tgPortfolio
rf_m = 0.01/12 mr_m = colMeans(return_m) mr.mtx_m = matrix(mr_m, ncol=1) mr_rf_m = mr_m - rf_m mr_rf_m = matrix(mr_rf_m, ncol=1) mr_rf_m library(fbasics) Sigma_m = cov(return_m) std_m = sqrt(return_m) ones = rep(1,3) one.vec = matrix(ones, ncol=1) a1_m = inv(Sigma_m)%%mr_rf_m b1_m = t(one.vec)%%a1_m tp_m = a1_m / as.numeric(b1_m) tp_m
ret.tp_m = sum(mr.mtx_m*tp_m) ret.tp_m
std.tp_m = sqrt((t(tp)%%Sigma_m)%%tp_m) std.tp_m
sharpe.tp_m = (ret.tp_m - rf_m)/std.tp_m sharpe.tp_m
ret.ts_w <- timeSeries(ret.w.dplyr,) chart.CumReturns(ret.ts_w, legend.loc = ‘topleft’, main = ’’)
plot(portfolioFrontier(ret.ts_w))
Spec = portfolioSpec() setSolver(Spec) = “solveRshortExact” setTargetReturn(Spec) = mean(colMeans(ret.ts_w))
Data_w <-ret.ts_w nAssets_w = getNAssets(portfolioData(Data_w)) Weights_w <- rep(1/nAssets_w, times = nAssets_w) covRisk(Data_w, Weights_w) varRisk(Data_w, Weights_w, alpha = 0.05)###VaR of equal weights portfolio cvarRisk(Data_w, Weights_w, alpha = 0.05)### CVaR for equal weight portfolio
return_w = ret.w.dplyr setRiskFreeRate(Spec) <- 0.01/52 tgPortfolio <- tangencyPortfolio(ret.ts_w, Spec, constraints = “Short”) tgPortfolio
rf_w = 0.01/52 mr_w = colMeans(return_w) mr.mtx_w = matrix(mr_w, ncol=1) mr_rf_w = mr_w - rf_w mr_rf_w = matrix(mr_rf_w, ncol=1) mr_rf_w library(fbasics) Sigma_w = cov(return_w) std_m = sqrt(return_m) ones = rep(1,3) one.vec = matrix(ones, ncol=1) a1_w = inv(Sigma_w)%%mr_rf_w b1_w = t(one.vec)%%a1_m tp_w = a1_m / as.numeric(b1_w) tp_w
ret.tp_w = sum(mr.mtx_w*tp_w) ret.tp_w
std.tp_w = sqrt((t(tp)%%Sigma_w)%%tp_w) std.tp_w
sharpe.tp_w = (ret.tp_w - rf_w)/std.tp_w sharpe.tp_w