Setsenchimeg
ch107035156
April 3, 2019
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.
# Reference:
# https://systematicinvestor.wordpress.com/category/factor-model/
# 3-asset example
#
rm(list=ls())
#setwd("D:/亞洲大學上課資料/Portfolio management 2017 Fall")
#
#=======================================
# Step 1: Import data from excel
#=======================================
# If data has percentage format, then you'd better convert it into numeric format in excel
# file = 3firmExample_data.csv
# First you have to convert factor data into numeric
#===========================================================================================================
# http://stackoverflow.com/questions/28977777/r-converting-data-frame-of-percentages-from-factor-to-numeric
#==========================================================================================================
# Delete % from data and convert into numeric
#==========================================================================================================
# After conversion, we save the file as 3firmExample_data3.csv
# setwd("~/portfolio_2017_Fall")
firm_data1 = read.csv('3firmExample_data3.csv')
str(firm_data1)## 'data.frame': 59 obs. of 4 variables:
## $ date : Factor w/ 59 levels "1995/10/1","1995/11/1",..: 4 5 6 7 8 9 10 1 2 3 ...
## $ Nordstrom: num -0.03615 -0.0568 0.07821 -0.00302 -0.02757 ...
## $ Starbucks: num 0.00521 -0.02105 0.21244 0.2036 0.04797 ...
## $ Microsoft: num 0.1213 0.13923 0.03529 0.06501 0.00138 ...firm_data1$date## [1] 1995/3/1 1995/4/1 1995/5/1 1995/6/1 1995/7/1 1995/8/1 1995/9/1
## [8] 1995/10/1 1995/11/1 1995/12/1 1996/1/1 1996/2/1 1996/3/1 1996/4/1
## [15] 1996/5/1 1996/6/1 1996/7/1 1996/8/1 1996/9/1 1996/10/1 1996/11/1
## [22] 1996/12/1 1997/1/1 1997/2/1 1997/3/1 1997/4/1 1997/5/1 1997/6/1
## [29] 1997/7/1 1997/8/1 1997/9/1 1997/10/1 1997/11/1 1997/12/1 1998/1/1
## [36] 1998/2/1 1998/3/1 1998/4/1 1998/5/1 1998/6/1 1998/7/1 1998/8/1
## [43] 1998/9/1 1998/10/1 1998/11/1 1998/12/1 1999/1/1 1999/2/1 1999/3/1
## [50] 1999/4/1 1999/5/1 1999/6/1 1999/7/1 1999/8/1 1999/9/1 1999/10/1
## [57] 1999/11/1 1999/12/1 2000/1/1
## 59 Levels: 1995/10/1 1995/11/1 1995/12/1 1995/3/1 1995/4/1 ... 2000/1/1library(xts)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(PerformanceAnalytics)##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legenddate1 = as.Date(firm_data1[,1], "%Y/%m/%d")
#convert firm_data1 into time series data: xts
firm_data1.xts = as.xts(firm_data1[,-1], order.by = date1)
firm.data1<-coredata(firm_data1.xts)
summary(firm.data1)## Nordstrom Starbucks Microsoft
## Min. :-0.212880 Min. :-0.47970 Min. :-0.17634
## 1st Qu.:-0.057395 1st Qu.:-0.01734 1st Qu.:-0.01826
## Median : 0.004950 Median : 0.04064 Median : 0.03848
## Mean : 0.001545 Mean : 0.02846 Mean : 0.04271
## 3rd Qu.: 0.064980 3rd Qu.: 0.09416 3rd Qu.: 0.11174
## Max. : 0.312480 Max. : 0.27967 Max. : 0.28153rbind(apply(firm.data1, 2, summary),
apply(firm.data1, 2, skewness),
apply(firm.data1, 2, kurtosis))## Nordstrom Starbucks Microsoft
## Min. -0.212880000 -0.47970000 -0.17634000
## 1st Qu. -0.057395000 -0.01734000 -0.01826500
## Median 0.004950000 0.04064000 0.03848000
## Mean 0.001545085 0.02846068 0.04271153
## 3rd Qu. 0.064980000 0.09415500 0.11173500
## Max. 0.312480000 0.27967000 0.28153000
## 0.242288454 -0.88822851 0.17120676
## 0.351952075 1.85144933 -0.08728437#====================================================================================
# IF you know the ticker of stocks, then you can
# download directly from yahoo finance
#=====================================================================================library(plyr)
library(quantmod)## Loading required package: TTR## Version 0.4-0 included new data defaults. See ?getSymbols.tickers<-c("JWN", "SBUX", "MSFT")
data.env<-new.env()
# here we use l_ply so that we don't double save the data
# getSymbols() does this already so we just want to be memory efficient
# go through every stock and try to use getSymbols()
l_ply(tickers, function(sym) try(getSymbols(sym, env=data.env), silent=T))## '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.# now we only want the stocks that got stored from getSymbols()
# basically we drop all "bad" tickers
stocks <- tickers[tickers %in% ls(data.env)]
# now we just loop through and merge our good stocks
# if you prefer to use an lapply version here, that is also fine
# since now we are just collecting all the good stock xts() objects
data <- xts()
# i=1
for(i in seq_along(stocks)) {
symbol <- stocks[i]
data <- merge(data, Ad(get(symbol, envir=data.env)))
}
head(data)## JWN.Adjusted SBUX.Adjusted MSFT.Adjusted
## 2007-01-03 08:00:00 35.71681 14.13297 22.47883
## 2007-01-04 08:00:00 36.32145 14.14900 22.44119
## 2007-01-05 08:00:00 35.64033 14.08886 22.31321
## 2007-01-08 08:00:00 35.82799 14.03674 22.53152
## 2007-01-09 08:00:00 36.21721 13.97660 22.55411
## 2007-01-10 08:00:00 36.39095 13.93249 22.32826str(data)## An 'xts' object on 2007-01-03 08:00:00/2019-04-02 08:00:00 containing:
## Data: num [1:3083, 1:3] 35.7 36.3 35.6 35.8 36.2 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:3] "JWN.Adjusted" "SBUX.Adjusted" "MSFT.Adjusted"
## Indexed by objects of class: [POSIXct,POSIXt] TZ:
## xts Attributes:
## NULL# convert POSIXct into date series
data<-xts(coredata(data), order.by = as.Date(index(data), tz=""))
head(data)## JWN.Adjusted SBUX.Adjusted MSFT.Adjusted
## 2007-01-03 35.71681 14.13297 22.47883
## 2007-01-04 36.32145 14.14900 22.44119
## 2007-01-05 35.64033 14.08886 22.31321
## 2007-01-08 35.82799 14.03674 22.53152
## 2007-01-09 36.21721 13.97660 22.55411
## 2007-01-10 36.39095 13.93249 22.32826library(fBasics)## Loading required package: timeDate##
## Attaching package: 'timeDate'## The following objects are masked from 'package:PerformanceAnalytics':
##
## kurtosis, skewness## Loading required package: timeSeries##
## Attaching package: 'timeSeries'## The following object is masked from 'package:zoo':
##
## time<-##
## Attaching package: 'fBasics'## The following object is masked from 'package:TTR':
##
## volatilitySigma = cov(firm_data1[,2:4])
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]
## Nordstrom 0.3635998
## Starbucks 0.1936537
## Microsoft 0.4427465mvp.ret<-sum((mvp.w)*colMeans(firm_data1[,2:4]))
mvp.ret## [1] 0.02498369#==================================
# Assume return is specified as 0.06.
# Try to find its optimal weight and standard deviation (tangency portfolio),
# expected return and Sharpe ratio.
#=================================mu<-0.06/12
return <- firm_data1[,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, 1)
out<-solve(Ax, b0)
wgt<-out[1:3]
wgt## Nordstrom Starbucks Microsoft
## 0.875635380 0.116816458 0.007548163sum(wgt)## [1] 1ret.out<-sum(wgt*colMeans(return))
ret.out.annual<-ret.out*12
ret.out.annual## [1] 0.06std.out<-sqrt(t(wgt)%*%cov(return)%*%wgt)
std.out.annual<-std.out*sqrt(12)
std.out.annual## [,1]
## [1,] 0.335302return = firm_data1[,2:4]
#specified portfolio return: mu
mu=0.06/12
minvariance <- function(return, mu) {
#return <- log(tail(assets, -1) / head(assets, -1))
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)
zx<-solve(Ax, b0)
weight<-zx[1:ncol(return)]
ret.out<-sum(weight*colMeans(return))
std.out<-sqrt(t(wgt)%*%cov(return)%*%wgt)
list(weight=weight, rtn=ret.out, sd=std.out)
}
minvariance(return, mu)## $weight
## Nordstrom Starbucks Microsoft
## 0.875635380 0.116816458 0.007548163
##
## $rtn
## [1] 0.005
##
## $sd
## [,1]
## [1,] 0.09679334frontier <- function(return){
#return <- log(tail(assets, -1) / head(assets, -1))
n = ncol(return)
Q = cov(return)
Ax <- rbind(2*cov(return), colMeans(return), rep(1, n))
Ax <- cbind(Ax, rbind(t(tail(Ax, 2)), matrix(0, 2, 2)))
r <- colMeans(return)
rbase <- seq(min(r), max(r), length = 100)
s <- sapply(rbase, function(x) {
b0 <- c(rep(0, ncol(return)), x, 1)
y <- head(solve(Ax, b0), n)
sqrt(y%*%Q%*%y)
})
plot(s, rbase, xlab = 'Std', ylab = 'Return')
}
frontier(return)
library(timeSeries)
library(PerformanceAnalytics)
#install.packages("rugarch", dependencies=TRUE)
#install.packages("PerformanceAnalytics", dependencies=TRUE)
#install.packages("fAssets", dependencies=TRUE)
#install.packages("fPortfolio",dependencies=TRUE)
return = firm_data1[,2:4]
# convert data to timeseries
ret.ts<- timeSeries(return, date1)
chart.CumReturns(ret.ts, legend.loc = 'topleft', main = '')
#To mimic what we have implemented in the preceding code, let us render the
#frontier plot of short sale constraints
#Spec = portfolioSpec()
#setSolver(Spec) = "solveRshortExact"
#setTargetReturn(Spec) = mean(colMeans(ret.ts))## or set your own target return
#DIFFERENT COVARIANCE ESTIMATORS
#1. MCd ESTIMATOR
#setEstimator(Spec)="covMcdEstimator"
#2. OGK estimator
#setEstimator(Spec)= "covOGKEstimator"
#3. Shrinkage estimator
#setEstimator(Spec)= "shrinkEstimator"
# Display structure of constraints;
#Spec
Constraints="Short"
#efficientPortfolio(ret.ts, Spec, Constraints)
#tangencyPortfolio(ret.ts, Spec, Constraints)
#
#portfolioConstraints(ret.ts, Spec, constraints = "LongOnly")
#setSolver(Spec) = "solveRquadprog"
#Frontier <- portfolioFrontier(as.timeSeries(ret.ts), Spec, constraints = "Short") #vs "LongOnly" "Short"
#frontierPlot(Frontier, col = c("orange", "red"), pch = 19)
#sharpeRatioLines(Frontier, col = "orange", lwd = 2) #### shows the Sharpe Ratio line
#monteCarloPoints(Frontier, mcSteps = 1000, cex = 0.25, pch = 19)
grid()
#
#weightsPlot(Frontier)
#EQUAL WEIGTH PORTFOLIO
Data<-ret.ts
#nAssets = getNAssets(portfolioData(Data))
#Weights <- rep(1/nAssets, times = nAssets)
#covRisk(Data, Weights)
#varRisk(Data, Weights, alpha = 0.05)###VaR of equal weights portfolio
#cvarRisk(Data, Weights, alpha = 0.05)### CVaR for equal weight portfoliorf = 0.01/12
mr = colMeans(firm_data1[,2:4])
mr.mtx = matrix(mr, ncol=1)
mr_rf = mr - rf
mr_rf = matrix(mr_rf, ncol=1)
mr_rf## [,1]
## [1,] 0.0007117514
## [2,] 0.0276273446
## [3,] 0.0418781921a1 = inv(Sigma)%*%mr_rf
b1 = t(one.vec)%*%a1
tp = a1 / as.numeric(b1)
tp## [,1]
## Nordstrom -0.2061492
## Starbucks 0.2791516
## Microsoft 0.9269976#portfolio expected return
ret.tp = sum(mr.mtx*tp)
ret.tp## [1] 0.04721981#portfolio standard deviation
std.tp = sqrt((t(tp)%*%Sigma)%*%tp)
std.tp## [,1]
## [1,] 0.1015897#sharpe ratio
sharpe.tp = (ret.tp - rf)/std.tp
sharpe.tp## [,1]
## [1,] 0.456606#=============================================
# Global min var portfolio
#=============================================
#globminSpec <- portfolioSpec()
# Answer will not change whether Constraints is "Short" or "LongOnly"
#globminPortfolio <- minvariancePortfolio(as.timeSeries(ret.ts), spec = #globminSpec, constraints = "LongOnly")
#print(globminPortfolio)
col = rampPalette(ncol(ret.ts), "purple2green")
#weights <- 100 * as.vector(getWeights(globminPortfolio))
names <- as.vector(names(ret.ts))
#barplot(height = weights, names.arg = names,
#horiz = TRUE, las = 1, col = col)
#title(main = "Weights of Global Min Variance Portfolio",
#xlab = "Weights %")
#=================================
# quadratic programming
#=================================
library(quadprog)
mu = apply(firm_data1[,2:4], 2, mean)
Amat = cbind(rep(1,3),mu) # set the constraints matrix
muP = seq(.01,.08,length=300) # set of 300 possible given returns
# for the expect portfolio return
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
#{
bvec = c(1,muP[i]) # constraint vector
result = solve.QP(Dmat=2*Sigma,dvec=rep(0,3),Amat=Amat,bvec=bvec,meq=2)
sdP[i] = sqrt(result$value)
weights[i,] = result$solution#postscript("3firmExample.ps",width=6,height=5)
pdf("3firmExample.pdf",width=6,height=5)
plot(sdP,muP,type="l",xlim=c(0,0.25),ylim=c(0,0.08),lty=1) # plot
# the efficient frontier (and inefficient frontier)
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 0 0lines(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],"Nordstrom",cex=1, pos=4)
text(std[2],mu[2],"Starbucks",cex=1, pos=4)
text(std[3],mu[3],"Microsoft",cex=1, pos=4)
graphics.off()Amat1 = cbind(rep(1,3),mu, diag(1,nrow=3)) # set the constraints matrix
t(Amat1)## Nordstrom Starbucks Microsoft
## 1.00000 1.0000 1.0000
## mu 0.00155 0.0285 0.0427
## 1.00000 0.0000 0.0000
## 0.00000 1.0000 0.0000
## 0.00000 0.0000 1.0000# 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.
muP1 = seq(min(mu)+.0001,max(mu)-.0001,length=300)
# for the expect portfolio return
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=2)
sdP1[i] = sqrt(result$value)
weights1[i,] = result$solution
}
#postscript("3firmExample_noshort.ps",width=6,height=5)
pdf("3firmExample_noshort.pdf",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## [1] 0 0 0lines(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],"Nordstrom",cex=1, pos=4)
text(std[2],mu[2],"Starbucks",cex=1, pos=4)
text(std[3],mu[3],"Microsoft",cex=1, pos=4)
graphics.off() create efficient frontier
ifelse(!require(corpcor), install.packages(“corpcor”), library(corpcor))
ifelse(!require(lpSolve), install.packages(“lpSolve”), library(lpSolve))
ef = portopt(ia, constraints, 50, ‘Efficient Frontier’)
ef
====================================================
David Ruppert’s example in his textbook
===================================================
ifelse(!require(Ecdat), install.packages(“Ecdat”), library(Ecdat))
ifelse(!require(Ecfun), install.packages(“Ecfun”), library(Ecfun))
ifelse(!require(quadprog), install.packages(“quadprog”), library(quadprog))
data(CRSPday) #daily observations from 1969-1-03 to 1998-12-31 #number of observations : 2528 #ge the return for General Electric, Permno 12060 #ibm the return for IBM, Permno 12490 #mobil the return for Mobil Corporation, Permno 15966 #crsp the return for the CRSP value-weighted index, including dividends class(CRSPday) CRSP.df = as.data.frame(CRSPday) head(CRSP.df) # R = 100*CRSP.df[,4:6] mean_vect = apply(R,2,mean) cov_mat = cov(R) sd_vect = sqrt(diag(cov_mat))
Amat = cbind(rep(1,3),mean_vect)
muP = seq(.05,.14,length=300) # set of 300 possible target values # for the expect portfolio return sdP = muP # set up storage for standard deviations of portfolio returns weights = matrix(0,nrow=300,ncol=3) # storage for portfolio weights
for (i in 1:length(muP)) # find the optimal portfolios for each target expected return { bvec = c(1,muP[i]) # constraint vector result = solve.QP(Dmat=2*cov_mat,dvec=rep(0,3),Amat=Amat,bvec=bvec,meq=2) sdP[i] = sqrt(result\(value) weights[i,] = result\)solution }
postscript(“CRSP_3firm.ps”,width=6,height=5) # Figure 11.3 plot(sdP,muP,type=“l”,xlim=c(0,2.5),ylim=c(0,.15),lty=3) # plot # the efficient frontier (and inefficient frontier) mufree = 1.3/253 # input value of risk-free interest rate points(0,mufree,cex=4,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=4,lty=2) # show line of optimal portfolios points(sdP[ind],muP[ind],cex=4,pch=“*“, col=‘blue’) # show tangency portfolio ind2 = (sdP == min(sdP)) # find the minimum variance portfolio points(sdP[ind2],muP[ind2],cex=2,pch=”+“) # show minimum variance portfolio ind3 = (muP > muP[ind2]) lines(sdP[ind3],muP[ind3],type=”l“,xlim=c(0,.25), ylim=c(0,.3),lwd=2) # plot the efficient frontier points(c(sd_vect[1],sd_vect[2], sd_vect[3]), c(mean_vect[1], mean_vect[2], mean_vect[3]), cex=1, pch=”o“, col=”red“) # show tangency portfolio text(sd_vect[1],mean_vect[1],”GE“,cex=1, pos=4) text(sd_vect[2],mean_vect[2],”IBM“,cex=1, pos=4) text(sd_vect[3],mean_vect[3],”Mobil“,cex=1, pos=4) graphics.off() ```