HW_01
Chantsal 108035487
4/3/2020
Q1. Try to find the min variance portfolio for three stocks
Using 2007-2009 daily returns as the insample data
library(quantmod)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: TTR## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Version 0.4-0 included new data defaults. See ?getSymbols.library(plyr)
tickers <- c("COST", "WMT", "TGT")
getSymbols(tickers, from= "2007-01-01", to= "2009-12-31")## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.## [1] "COST" "WMT" "TGT"data.env<-new.env()
l_ply(tickers, function(sym) try(getSymbols(sym, env=data.env), silent=T))
stocks <- tickers[tickers %in% ls(data.env)]
data <- xts()
for(i in seq_along(stocks)) {
symbol <- stocks[i]
data <- merge(data, Ad(get(symbol, envir=data.env)))
}
head(data)## COST.Adjusted WMT.Adjusted TGT.Adjusted
## 2007-01-03 39.56923 34.86750 41.24798
## 2007-01-04 40.49031 35.03615 41.54374
## 2007-01-05 39.99606 34.75017 41.31289
## 2007-01-08 40.16081 34.46419 41.37062
## 2007-01-09 40.46037 34.75017 42.02705
## 2007-01-10 40.74491 34.66951 42.06313str(data)## An 'xts' object on 2007-01-03/2020-04-03 containing:
## Data: num [1:3337, 1:3] 39.6 40.5 40 40.2 40.5 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:3] "COST.Adjusted" "WMT.Adjusted" "TGT.Adjusted"
## Indexed by objects of class: [POSIXct,POSIXt] TZ: Etc/UTC
## xts Attributes:
## NULLcon = gzcon(url('https://github.com/systematicinvestor/SIT/raw/master/sit.gz', 'rb'))
source(con)
close(con)In other words, you have to compute optimal weights for 2010-01
Also given the weights from Q1, compute realized returns for 2010-01.library(fBasics)## Loading required package: timeDate## Loading required package: timeSeries##
## Attaching package: 'timeSeries'## The following object is masked from 'package:zoo':
##
## time<-##
## Attaching package: 'fBasics'## The following objects are masked _by_ '.GlobalEnv':
##
## inv, vec## The following object is masked from 'package:TTR':
##
## volatilitylibrary(xts)
library(PerformanceAnalytics)##
## Attaching package: 'PerformanceAnalytics'## The following objects are masked from 'package:timeDate':
##
## kurtosis, skewness## The following object is masked from 'package:graphics':
##
## legendfirm.data= data
str(firm.data)## An 'xts' object on 2007-01-03/2020-04-03 containing:
## Data: num [1:3337, 1:3] 39.6 40.5 40 40.2 40.5 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:3] "COST.Adjusted" "WMT.Adjusted" "TGT.Adjusted"
## Indexed by objects of class: [POSIXct,POSIXt] TZ: Etc/UTC
## xts Attributes:
## NULLSigma = cov(firm.data)
std = sqrt(diag(Sigma))
ones = rep(1,3)
one.vec = matrix(ones, ncol=1)
a = inv(Sigma)%*%one.vec
b = t(one.vec)%*%a
mvp.w =a / as.numeric(b)
mvp.w## [,1]
## COST.Adjusted -0.3681164
## WMT.Adjusted 0.7702422
## TGT.Adjusted 0.5978741mvp.ret<-sum((mvp.w)*colMeans(firm.data))
mvp.ret## [1] 37.98572Give me weights, realized returns and standard deviation for 2010-01.
mu<-0.05/12
return <- firm.data
Ax <- rbind(2*cov(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## COST.Adjusted WMT.Adjusted TGT.Adjusted
## -0.90108684 -0.09534515 1.99643199sum(wgt)## [1] 1ret.out<-sum(wgt*colMeans(return))
ret.out.annual<-ret.out*12
ret.out.annual## [1] 0.05std.out<-sqrt(t(wgt)%*%cov(return)%*%wgt)
std.out.annual<-std.out*sqrt(12)
std.out.annual## [,1]
## [1,] 119.8548Q2. Based on Q1, you have to use weekly and monthly returns to get the answers to Q1
rm(list=ls())
con = gzcon(url('https://github.com/systematicinvestor/SIT/raw/master/sit.gz', 'rb'))
source(con)
close(con)
library(quantmod)
tickers <- c("COST","WMT","TGT")
n <- length(tickers)
data <- new.env()
getSymbols(tickers, src = 'yahoo', from = '2007-01-01',to= '2009-12-31', env = data, auto.assign = T)## [1] "COST" "WMT" "TGT"for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)
data.monthly <- new.env()
for(i in tickers) data.monthly[[i]] = to.monthly(data[[i]], indexAt='endof')
data.monthly## <environment: 0x555a3af499c8>bt.prep(data, align='remove.na', fill.gaps = T)
names(data)## [1] "prices" "COST" "dates" "WMT"
## [5] "weight" ".getSymbols" "TGT" "symbolnames"
## [9] "execution.price"head(data$prices)## COST TGT WMT
## 2007-01-03 39.56923 41.24798 34.86750
## 2007-01-04 40.49031 41.54374 35.03615
## 2007-01-05 39.99606 41.31289 34.75017
## 2007-01-08 40.16081 41.37062 34.46419
## 2007-01-09 40.46037 42.02705 34.75017
## 2007-01-10 40.74491 42.06313 34.66951head(data$weight)## COST TGT WMT
## 2007-01-03 NA NA NA
## 2007-01-04 NA NA NA
## 2007-01-05 NA NA NA
## 2007-01-08 NA NA NA
## 2007-01-09 NA NA NA
## 2007-01-10 NA NA NAprices <- data$prices
ret <- prices/mlag(prices)- 1
head(ret)## COST TGT WMT
## 2007-01-03 NA NA NA
## 2007-01-04 0.023277631 0.0071702174 0.004837055
## 2007-01-05 -0.012206747 -0.0055568903 -0.008162340
## 2007-01-08 0.004119306 0.0013974333 -0.008229800
## 2007-01-09 0.007458838 0.0158672277 0.008298092
## 2007-01-10 0.007032809 0.0008584708 -0.002321398ret<- ret[-1, ]
head(ret)## COST TGT WMT
## 2007-01-04 0.023277631 0.0071702174 0.004837055
## 2007-01-05 -0.012206747 -0.0055568903 -0.008162340
## 2007-01-08 0.004119306 0.0013974333 -0.008229800
## 2007-01-09 0.007458838 0.0158672277 0.008298092
## 2007-01-10 0.007032809 0.0008584708 -0.002321398
## 2007-01-11 0.016725056 0.0257240705 0.006768311Use weekly
library(tibbletime)##
## Attaching package: 'tibbletime'## The following object is masked from 'package:timeSeries':
##
## filter## The following object is masked from 'package:stats':
##
## filterlibrary(timetk)## Loading required package: recipes## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked _by_ '.GlobalEnv':
##
## count, lst## The following objects are masked from 'package:timeSeries':
##
## filter, lag## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following objects are masked from 'package:xts':
##
## first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union##
## Attaching package: 'recipes'## The following object is masked from 'package:stats':
##
## stepweekly <- prices %>%
tk_tbl(rename_index= "date") %>%
as_tbl_time(index= date) %>%
as_period(period= "week", side= "end")
#Q1. min variance
str(weekly[,2:4])## Classes 'tbl_df', 'tbl' and 'data.frame': 157 obs. of 3 variables:
## $ COST: num 40 41.9 42.5 41.5 41.9 ...
## $ TGT : num 41.3 43.5 45 43.4 44.7 ...
## $ WMT : num 34.8 35.2 35.4 35 35.3 ...Sigma.w = cov(weekly[,2:4])
std = sqrt(diag(Sigma.w))
ones = rep(1,3)
one.vec = matrix(ones, ncol=1)
a = inv(Sigma.w)%*%one.vec
b = t(one.vec)%*%a
mvp.w.w =a / as.numeric(b)
mvp.w.w## [,1]
## COST -0.2188173
## TGT 0.3645299
## WMT 0.8542874mvp.ret.w<-sum((mvp.w.w)*colMeans(weekly[,2:4]))
mvp.ret.w## [1] 36.39087mu.w<-0.05/12
return<- weekly[,2:4]
Ax <- rbind(2*cov(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.w, 1)
out.w<-solve(Ax, b0)
wgt<-out.w[1:3]
wgt## COST TGT WMT
## -5.565384 3.195241 3.370143sum(wgt)## [1] 1ret.out<-sum(wgt*colMeans(return))
ret.out.annual<-ret.out*12
ret.out.annual## [1] 0.05std.out<-sqrt(t(wgt)%*%cov(return)%*%wgt)
std.out.annual<-std.out*sqrt(12)
std.out.annual## [,1]
## [1,] 81.57326use montly data
monthly <- prices %>%
tk_tbl(rename_index= "date") %>%
as_tbl_time(index= date) %>%
as_period(period= "month", side= "end")
#Q1. min variance
str(monthly[,2:4])## Classes 'tbl_df', 'tbl' and 'data.frame': 36 obs. of 3 variables:
## $ COST: num 42.1 42 40.4 40.3 42.5 ...
## $ TGT : num 44.3 44.5 42.8 42.9 45.2 ...
## $ WMT : num 35 35.4 34.6 35.3 35.2 ...Sigma.m = cov(monthly[,2:4])
std = sqrt(diag(Sigma.m))
ones = rep(1,3)
one.vec = matrix(ones, ncol=1)
a = inv(Sigma.m)%*%one.vec
b = t(one.vec)%*%a
mvp.w.m =a / as.numeric(b)
mvp.w.m## [,1]
## COST -0.2573858
## TGT 0.4122964
## WMT 0.8450894mvp.ret.m<-sum((mvp.w.m)*colMeans(monthly[,2:4]))
mvp.ret.m## [1] 36.16847mu.m<-0.05/12
return.m<-monthly[,2:4]
Ax.m <- rbind(2*cov(return.m), colMeans(return.m), rep(1, ncol(return.m)))
Ax.m <- cbind(Ax, rbind(t(tail(Ax.m, 2)), matrix(0, 2, 2)))
b0 <- c(rep(0, ncol(return.m)), mu.m, 1)
out.m<-solve(Ax, b0)
wgt.m<-out.m[1:3]
wgt.m## COST TGT WMT
## -5.565384 3.195241 3.370143sum(wgt.m)## [1] 1ret.out.m<-sum(wgt.m*colMeans(return.m))
ret.out.annual.m<-ret.out.m*12
ret.out.annual.m## [1] -1.700172std.out<-sqrt(t(wgt)%*%cov(return)%*%wgt)
std.out.annual<-std.out*sqrt(12)
std.out.annual## [,1]
## [1,] 81.57326Q3. Find the tangency portfolio based on Q2.
mean <- apply(ret, 2, mean)
mean## COST TGT WMT
## 0.0004121230 0.0002460882 0.0003854028var <- apply(ret, 2, var)
var## COST TGT WMT
## 0.0003964414 0.0008174907 0.0002671032summary <- apply(ret, 2, summary)
apply(ret, 2, skewness)## COST TGT WMT
## 0.2087893 0.2928760 0.4343346cov.matrix <- matrix(cov(ret), ncol =3)
cov.matrix## [,1] [,2] [,3]
## [1,] 0.0003964414 0.0003702955 0.0002106663
## [2,] 0.0003702955 0.0008174907 0.0002875378
## [3,] 0.0002106663 0.0002875378 0.0002671032mu <- matrix(mean, ncol=1)
mu## [,1]
## [1,] 0.0004121230
## [2,] 0.0002460882
## [3,] 0.0003854028rf <- 0.005
numerator <- solve(cov.matrix)%*% (mu - rf*one.vec)
denominator <-t(one.vec)%*% numerator
tangency.weight <-numerator/as.numeric(denominator)
tangency.return<- t(mu)%*% tangency.weight
tangency.var <- t(tangency.weight) %*% cov.matrix %*%tangency.weight
tangency.std <- sqrt(as.numeric(tangency.var))
Sharperatio <- (tangency.return/tangency.std)
Sharperatio## [,1]
## [1,] 0.02589659#HW_02 ## Q1. Find the tangency portfolio with no short-sale constraint. Risk-free rate is assumed to be 1% per annum.
library(quantmod)
library(plyr)
library(xts)
library(quadprog)
library(fBasics)
library(stats)
tickers <- c("COST", "WMT", "TGT")
getSymbols(tickers, from= "2007-01-01", to= "2009-12-31")## [1] "COST" "WMT" "TGT"data.env<-new.env()
l_ply(tickers, function(sym) try(getSymbols(sym, env=data.env), silent=T))
stocks <- tickers[tickers %in% ls(data.env)]
data <- xts()
for(i in seq_along(stocks)) {
symbol <- stocks[i]
data <- merge(data, Ad(get(symbol, envir=data.env)))
}
mu = apply(data, 2, mean)
Sigma = cov(data)
std = sqrt(diag(Sigma))
Amat1 = cbind(rep(1,3),mu, diag(1,nrow=3))
t(Amat1)## COST.Adjusted WMT.Adjusted TGT.Adjusted
## 1.0000 1.00000 1.00000
## mu 113.1966 61.47748 54.02924
## 1.0000 0.00000 0.00000
## 0.0000 1.00000 0.00000
## 0.0000 0.00000 1.00000# muP = seq(.01,.08,length=300) # set of 300 possible target values
# When short sales are prohibited, the target expected return on the
# portfolio must lie between the smallest
# and largest expected returns on the stocks.
muP = seq(min(mu)+.0001,max(mu)-.0001,length=300)
sdP = muP # set up storage for standard deviations of portfolio returns
weights = matrix(0,nrow=300,ncol=3) # storage for portfolio weights
i=1
for (i in 1:length(muP)) # find the optimal portfolios for each target expected return
{
bvec1 = c(1,muP[i], rep(0,3)) # constraint vector
result = solve.QP(Dmat=2*Sigma,dvec=rep(0,3),Amat=Amat1,bvec=bvec1,meq=2)
sdP[i] = sqrt(result$value)
weights[i,] = result$solution
}
plot(sdP,muP,type="l",xlim=c(0,0.25),ylim=c(0,0.08),lty=1)
par(new=TRUE)
plot(sdP,muP,type="l",xlim=c(0,0.25),ylim=c(0,0.08),lty=1, col="green")
mufree = 0.01/12 # 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## [1] 0.000 0.372 0.628lines(c(0,2),mufree+c(0,2)*(muP[ind]-mufree)/sdP[ind],lwd=2,lty=3)
# show line of optimal portfolios
points(sdP[ind],muP[ind],cex=1,pch=19, col="red") # show tangency portfolio
ind2 = (sdP == min(sdP)) # find the minimum variance portfolio
points(sdP[ind2],muP[ind2],cex=2,pch="+", col="blue") # show minimum variance portfolio
ind3 = (muP > muP[ind2])
lines(sdP[ind3],muP[ind3],type="l",xlim=c(0,.25),
ylim=c(0,.08),lwd=1) # plot the efficient frontier
points(c(std[1],std[2], std[3]), c(mu[1], mu[2], mu[3]), cex=1, pch="o", col="red")
text(std[1],mu[1],"Costco",cex=1, pos=4)
text(std[2],mu[2],"Wallmart",cex=1, pos=4)
text(std[3],mu[3],"Target",cex=1, pos=4)
graphics.off()Q2.Plot the efficient frontier with and without short-sale constraint.
con = gzcon(url('https://github.com/systematicinvestor/SIT/raw/master/sit.gz', 'rb'))
source(con)
close(con)
tickers <- c("COST", "WMT", "TGT")
n <- length(tickers)
#with short sale constraint.
constraints = new.constraints(n, lb = -Inf, ub = +Inf)
constraints = add.constraints(rep(1, n), 1, type = '=', constraints)
ia<- create.historical.ia(ret, 250)
weight <- min.risk.portfolio(ia, constraints)
#without short sale constraint.
ia<- create.historical.ia(ret, 250)
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)
weight## [1] 2.33e-01 1.39e-17 7.67e-01ifelse(!require(corpcor), install.packages("corpcor"), library(corpcor))## Loading required package: corpcor##
## Attaching package: 'corpcor'## The following object is masked _by_ '.GlobalEnv':
##
## cov.shrink## [1] "corpcor"ifelse(!require(lpSolve), install.packages("lpSolve"), library(lpSolve))## Loading required package: lpSolve## [1] "lpSolve"ef = portopt(ia, constraints, 50, 'Efficient Frontier') ##
## Attaching package: 'kernlab'## The following object is masked _by_ '.GlobalEnv':
##
## crossef## $weight
## COST TGT WMT
## [1,] 0.233 1.39e-17 0.7670
## [2,] 0.249 -1.30e-17 0.7513
## [3,] 0.264 -9.57e-18 0.7357
## [4,] 0.280 -6.16e-18 0.7200
## [5,] 0.296 -2.75e-18 0.7044
## [6,] 0.311 6.66e-19 0.6887
## [7,] 0.327 4.08e-18 0.6731
## [8,] 0.343 7.49e-18 0.6574
## [9,] 0.358 1.09e-17 0.6418
## [10,] 0.374 1.43e-17 0.6261
## [11,] 0.390 1.77e-17 0.6105
## [12,] 0.405 2.11e-17 0.5948
## [13,] 0.421 2.45e-17 0.5792
## [14,] 0.436 2.80e-17 0.5635
## [15,] 0.452 3.14e-17 0.5479
## [16,] 0.468 3.48e-17 0.5322
## [17,] 0.483 3.82e-17 0.5165
## [18,] 0.499 4.16e-17 0.5009
## [19,] 0.515 4.50e-17 0.4852
## [20,] 0.530 4.84e-17 0.4696
## [21,] 0.546 5.18e-17 0.4539
## [22,] 0.562 5.53e-17 0.4383
## [23,] 0.577 5.87e-17 0.4226
## [24,] 0.593 6.21e-17 0.4070
## [25,] 0.609 6.55e-17 0.3913
## [26,] 0.624 6.89e-17 0.3757
## [27,] 0.640 7.23e-17 0.3600
## [28,] 0.656 7.57e-17 0.3444
## [29,] 0.671 7.91e-17 0.3287
## [30,] 0.687 8.26e-17 0.3131
## [31,] 0.703 8.60e-17 0.2974
## [32,] 0.718 8.94e-17 0.2818
## [33,] 0.734 9.28e-17 0.2661
## [34,] 0.750 9.62e-17 0.2504
## [35,] 0.765 9.96e-17 0.2348
## [36,] 0.781 1.03e-16 0.2191
## [37,] 0.797 1.06e-16 0.2035
## [38,] 0.812 1.10e-16 0.1878
## [39,] 0.828 1.13e-16 0.1722
## [40,] 0.843 1.17e-16 0.1565
## [41,] 0.859 1.20e-16 0.1409
## [42,] 0.875 1.23e-16 0.1252
## [43,] 0.890 1.27e-16 0.1096
## [44,] 0.906 1.30e-16 0.0939
## [45,] 0.922 1.34e-16 0.0783
## [46,] 0.937 1.37e-16 0.0626
## [47,] 0.953 1.41e-16 0.0470
## [48,] 0.969 1.44e-16 0.0313
## [49,] 0.984 1.47e-16 0.0157
## [50,] 1.000 0.00e+00 0.0000
##
## $return
## [,1]
## [1,] 0.103
## [2,] 0.103
## [3,] 0.103
## [4,] 0.103
## [5,] 0.103
## [6,] 0.103
## [7,] 0.104
## [8,] 0.104
## [9,] 0.104
## [10,] 0.104
## [11,] 0.104
## [12,] 0.104
## [13,] 0.104
## [14,] 0.104
## [15,] 0.104
## [16,] 0.105
## [17,] 0.105
## [18,] 0.105
## [19,] 0.105
## [20,] 0.105
## [21,] 0.105
## [22,] 0.105
## [23,] 0.105
## [24,] 0.105
## [25,] 0.106
## [26,] 0.106
## [27,] 0.106
## [28,] 0.106
## [29,] 0.106
## [30,] 0.106
## [31,] 0.106
## [32,] 0.106
## [33,] 0.107
## [34,] 0.107
## [35,] 0.107
## [36,] 0.107
## [37,] 0.107
## [38,] 0.107
## [39,] 0.107
## [40,] 0.107
## [41,] 0.107
## [42,] 0.108
## [43,] 0.108
## [44,] 0.108
## [45,] 0.108
## [46,] 0.108
## [47,] 0.108
## [48,] 0.108
## [49,] 0.108
## [50,] 0.109
##
## $risk
## [1] 0.252 0.252 0.252 0.252 0.252 0.253 0.253 0.253 0.254 0.254 0.255 0.256
## [13] 0.256 0.257 0.258 0.259 0.259 0.260 0.261 0.262 0.263 0.265 0.266 0.267
## [25] 0.268 0.270 0.271 0.273 0.274 0.276 0.277 0.279 0.281 0.282 0.284 0.286
## [37] 0.288 0.289 0.291 0.293 0.295 0.297 0.299 0.302 0.304 0.306 0.308 0.310
## [49] 0.313 0.315
##
## $name
## [1] "Efficient Frontier"Q3. Repeat Q1 and Q2 using weekly and monthly returns.
library(tibbletime)
library(timetk)
library(fBasics)
library(quadprog)
weekly <- prices %>%
tk_tbl(rename_index= "date") %>%
as_tbl_time(index= date) %>%
as_period(period= "week", side= "end")
Sigma1 = cov(weekly[,2:4])
std1 = sqrt(diag(Sigma1))
mu = apply(weekly[,2:4], 2, mean)
Amat2 = cbind(rep(1,3),mu, diag(1,nrow=3))
t(Amat2)## COST TGT WMT
## 1.0 1.0 1.0
## mu 44.2 36.5 38.3
## 1.0 0.0 0.0
## 0.0 1.0 0.0
## 0.0 0.0 1.0muP2 = seq(min(mu)+.0001, max(mu)-.0001, length = 150)
sdP2 = muP2
weights2 = matrix(0,nrow=150,ncol=3)
i=1
for (i in 1:length(muP2))
{
bvec2 = c(1,muP2[i], rep(0,3))
result = solve.QP(Dmat=2*Sigma1,dvec=rep(0,3),Amat=Amat2,bvec=bvec2,meq=2)
sdP2[i] = sqrt(result$value)
weights2[i,] = result$solution
}
plot(sdP2,muP2,type="l",xlim=c(0,0.25),ylim=c(0,0.08),lty=1)
mufree.w = 0.01/12
points(0,mufree.w,cex=3, pch="*")
sharpe =( muP2-mufree.w)/sdP2
ind1 = (sharpe == max(sharpe))
options(digits=3)
weights2[ind1,] ## [1] 0.000482 0.248423 0.751095lines(c(0,2),mufree.w+c(0,2)*(muP2[ind1]-mufree.w)/sdP2[ind1],lwd=2,lty=3)
points(sdP2[ind1],muP2[ind1],cex=1,pch=19, col="red")
ind2 = (sdP2 == min(sdP2))
points(sdP2[ind2],muP2[ind2],cex=2,pch="+", col="blue")
ind3 = (muP2 > muP2[ind2])
lines(sdP2[ind3],muP2[ind3],type="l",xlim=c(0,.25),
ylim=c(0,.08),lwd=1)
points(c(std1[1],std1[2], std1[3]), c(mu[1], mu[2], mu[3]), cex=1, pch="o", col="red")
text(std1[1],mu[1],"Costco",cex=1, pos=4)
text(std1[2],mu[2],"Wallmart",cex=1, pos=4)
text(std1[3],mu[3],"Target",cex=1, pos=4)
graphics.off()monthly <- prices %>%
tk_tbl(rename_index= "date") %>%
as_tbl_time(index= date) %>%
as_period(period= "month", side= "end")
Sigma2 = cov(monthly[,2:4])
std2= sqrt(diag(Sigma2))
mu = apply(monthly[,2:4], 2, mean)
Amat3 = cbind(rep(1,3),mu, diag(1,nrow=3))
t(Amat3)## COST TGT WMT
## 1.0 1.0 1.0
## mu 44.3 36.5 38.5
## 1.0 0.0 0.0
## 0.0 1.0 0.0
## 0.0 0.0 1.0muP3 = seq(min(mu)+.0001,max(mu)-.0001,length=300)
sdP3= muP3
weights3= matrix(0,nrow=300,ncol=3)
i=1
for (i in 1:length(muP3))
{
bvec3 = c(1,muP3[i], rep(0,3))
result = solve.QP(Dmat=2*Sigma2,dvec=rep(0,3),Amat=Amat3,bvec=bvec3,meq=2)
sdP3[i] = sqrt(result$value)
weights3[i,] = result$solution
}
plot(sdP3,muP3,type="l",xlim=c(0,0.25),ylim=c(0,0.08),lty=1)
mufree.m = 0.01/12
points(0,mufree.m,cex=3, pch="*")
sharpe =( muP3-mufree.m)/sdP3
ind2 = (sharpe == max(sharpe))
options(digits=3)
weights3[ind2,] ## [1] -1.3e-18 2.8e-01 7.2e-01lines(c(0,2),mufree.m+c(0,2)*(muP3[ind2]-mufree.m)/sdP3[ind2],lwd=2,lty=3)
points(sdP3[ind2],muP3[ind2],cex=1,pch=19, col="red")
ind2 = (sdP3 == min(sdP3))
points(sdP3[ind2],muP3[ind2],cex=2,pch="+", col="blue")
ind3 = (muP3 > muP3[ind2])
lines(sdP3[ind3],muP3[ind3],type="l",xlim=c(0,.25),
ylim=c(0,.08),lwd=1)
points(c(std[1],std[2], std[3]), c(mu[1], mu[2], mu[3]), cex=1, pch="o", col="red")
text(std[1],mu[1],"Costco",cex=1, pos=4)
text(std[2],mu[2],"Wallmart",cex=1, pos=4)
text(std[3],mu[3],"Target",cex=1, pos=4)
graphics.off()