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.06313
str(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:
## NULL
con = 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':
##
## volatility
library(xts)
library(PerformanceAnalytics)
##
## Attaching package: 'PerformanceAnalytics'
## The following objects are masked from 'package:timeDate':
##
## kurtosis, skewness
## The following object is masked from 'package:graphics':
##
## legend
firm.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:
## NULL
Sigma = 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.5978741
mvp.ret<-sum((mvp.w)*colMeans(firm.data))
mvp.ret
## [1] 37.98572
Give 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.99643199
sum(wgt)
## [1] 1
ret.out<-sum(wgt*colMeans(return))
ret.out.annual<-ret.out*12
ret.out.annual
## [1] 0.05
std.out<-sqrt(t(wgt)%*%cov(return)%*%wgt)
std.out.annual<-std.out*sqrt(12)
std.out.annual
## [,1]
## [1,] 119.8548
Q2. 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: 0x56343aea5010>
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.66951
head(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 NA
prices <- 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.002321398
ret<- 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.006768311
Use weekly
library(tibbletime)
##
## Attaching package: 'tibbletime'
## The following object is masked from 'package:timeSeries':
##
## filter
## The following object is masked from 'package:stats':
##
## filter
library(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':
##
## step
weekly <- 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.8542874
mvp.ret.w<-sum((mvp.w.w)*colMeans(weekly[,2:4]))
mvp.ret.w
## [1] 36.39087
mu.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.370143
sum(wgt)
## [1] 1
ret.out<-sum(wgt*colMeans(return))
ret.out.annual<-ret.out*12
ret.out.annual
## [1] 0.05
std.out<-sqrt(t(wgt)%*%cov(return)%*%wgt)
std.out.annual<-std.out*sqrt(12)
std.out.annual
## [,1]
## [1,] 81.57326
use 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.8450894
mvp.ret.m<-sum((mvp.w.m)*colMeans(monthly[,2:4]))
mvp.ret.m
## [1] 36.16847
mu.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.370143
sum(wgt.m)
## [1] 1
ret.out.m<-sum(wgt.m*colMeans(return.m))
ret.out.annual.m<-ret.out.m*12
ret.out.annual.m
## [1] -1.700172
std.out<-sqrt(t(wgt)%*%cov(return)%*%wgt)
std.out.annual<-std.out*sqrt(12)
std.out.annual
## [,1]
## [1,] 81.57326
Q3. Find the tangency portfolio based on Q2.
mean <- apply(ret, 2, mean)
mean
## COST TGT WMT
## 0.0004121230 0.0002460882 0.0003854028
var <- apply(ret, 2, var)
var
## COST TGT WMT
## 0.0003964414 0.0008174907 0.0002671032
summary <- apply(ret, 2, summary)
apply(ret, 2, skewness)
## COST TGT WMT
## 0.2087893 0.2928760 0.4343346
cov.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.0002671032
mu <- matrix(mean, ncol=1)
mu
## [,1]
## [1,] 0.0004121230
## [2,] 0.0002460882
## [3,] 0.0003854028
rf <- 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