rm(list=ls())
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/1
library(xts)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(PerformanceAnalytics)
##
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
##
## legend
date1 = 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.28153
skewness(firm.data1)
## Nordstrom Starbucks Microsoft
## Skewness 0.2422885 -0.8882285 0.1712068
rbind(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(quantmod)
#tickers<-c("JWN", "SBUX", "MSFT")
#getSymbols(tickers, from = '2010-12-31', to = '2016-12-31', auto.assign = TRUE)
#=================================
# Minimum variance portfolio
#=================================
library(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<-
Sigma = 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.4427465
mvp.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.007548163
sum(wgt)
## [1] 1
ret.out<-sum(wgt*colMeans(return))
ret.out.annual<-ret.out*12
ret.out.annual
## [1] 0.06
std.out<-sqrt(t(wgt)%*%cov(return)%*%wgt)
std.out.annual<-std.out*sqrt(12)
std.out.annual
## [,1]
## [1,] 0.335302
firm_data1 = read.csv('3firmExample_data3.csv')
#==================================
# 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.007548163
sum(wgt)
## [1] 1
ret.out<-sum(wgt*colMeans(return))
ret.out.annual<-ret.out*12
ret.out.annual
## [1] 0.06
std.out<-sqrt(t(wgt)%*%cov(return)%*%wgt)
std.out.annual<-std.out*sqrt(12)
std.out.annual
## [,1]
## [1,] 0.335302
#====================================================================================================
# 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 = 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.09679334
#======================================================
# Create frontier function to plot efficient frontier
#======================================================
frontier <- 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)
#===============================================
# Use package fPortfolio to plot frontier
#===============================================
library(timeSeries)
library(PerformanceAnalytics)
#install.packages("rugarch", dependencies=TRUE)
#install.packages("PerformanceAnalytics", dependencies=TRUE)
#install.packages("fAssets", dependencies=TRUE)
#install.packages("fPortfolio",dependencies=TRUE)
library(fPortfolio)
## Loading required package: fAssets
return = firm_data1[,2:4]
date1 = as.Date(firm_data1[,1], "%Y/%m/%d")
ret.ts<- timeSeries(return, date1)
chart.CumReturns(ret.ts, legend.loc = 'topleft', main = '')
1#To mimic what we have implemented in the preceding code, let us render the
## [1] 1
#frontier plot of short sale constraints
Spec = portfolioSpec()
setSolver(Spec) = "solveRshortExact"
#setTargetReturn(Spec) = mean(colMeans(ret.ts))## or set your own target return
Spec
##
## Model List:
## Type: MV
## Optimize: minRisk
## Estimator: covEstimator
## Params: alpha = 0.05
##
## Portfolio List:
## Target Weights: NULL
## Target Return: NULL
## Target Risk: NULL
## Risk-Free Rate: 0
## Number of Frontier Points: 50
##
## Optim List:
## Solver: solveRshortExact
## Objective: portfolioObjective portfolioReturn portfolioRisk
## Trace: FALSE
#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
##
## Model List:
## Type: MV
## Optimize: minRisk
## Estimator: shrinkEstimator
## Params: alpha = 0.05
##
## Portfolio List:
## Target Weights: NULL
## Target Return: NULL
## Target Risk: NULL
## Risk-Free Rate: 0
## Number of Frontier Points: 50
##
## Optim List:
## Solver: solveRshortExact
## Objective: portfolioObjective portfolioReturn portfolioRisk
## Trace: FALSE
Constraints="Short"
efficientPortfolio(ret.ts, Spec, Constraints)
## Warning in as.vector(invSigma %*% one)/(one %*% invSigma %*% one): Recycling array of length 1 in vector-array arithmetic is deprecated.
## Use c() or as.vector() instead.
##
## Title:
## MV Efficient Portfolio
## Estimator: shrinkEstimator
## Solver: solveRshortExact
## Optimize: minRisk
## Constraints: Short
##
## Portfolio Weights:
## Nordstrom Starbucks Microsoft
## 0.3750 0.2047 0.4203
##
## Covariance Risk Budgets:
## Nordstrom Starbucks Microsoft
## 0.3823 0.2131 0.4045
##
## Target Returns and Risks:
## mean mu Cov Sigma CVaR VaR
## 0.0244 0.0244 0.0734 0.0667 0.1324 0.1237
##
## Description:
## Fri Mar 29 19:11:50 2019 by user:
tangencyPortfolio(ret.ts, Spec, Constraints)
##
## Title:
## MV Tangency Portfolio
## Estimator: shrinkEstimator
## Solver: solveRshortExact
## Optimize: minRisk
## Constraints: Short
##
## Portfolio Weights:
## Nordstrom Starbucks Microsoft
## -0.0147 0.2486 0.7661
##
## Covariance Risk Budgets:
## Nordstrom Starbucks Microsoft
## -0.0037 0.1924 0.8113
##
## Target Returns and Risks:
## mean mu Cov Sigma CVaR VaR
## 0.0398 0.0398 0.0870 0.0853 0.1485 0.1131
##
## Description:
## Fri Mar 29 19:11:50 2019 by user:
minvariancePortfolio(ret.ts, Spec, Constraints)
## Warning in as.vector(invSigma %*% one)/(one %*% invSigma %*% one): Recycling array of length 1 in vector-array arithmetic is deprecated.
## Use c() or as.vector() instead.
##
## Title:
## MV Minimum Variance Portfolio
## Estimator: shrinkEstimator
## Solver: solveRshortExact
## Optimize: minRisk
## Constraints: Short
##
## Portfolio Weights:
## Nordstrom Starbucks Microsoft
## 0.3750 0.2047 0.4203
##
## Covariance Risk Budgets:
## Nordstrom Starbucks Microsoft
## 0.3823 0.2131 0.4045
##
## Target Returns and Risks:
## mean mu Cov Sigma CVaR VaR
## 0.0244 0.0244 0.0734 0.0667 0.1324 0.1237
##
## Description:
## Fri Mar 29 19:11:50 2019 by user:
#
portfolioConstraints(ret.ts, Spec, constraints = "LongOnly")
##
## Title:
## Portfolio Constraints
##
## Lower/Upper Bounds:
## Nordstrom Starbucks Microsoft
## Lower 0 0 0
## Upper 1 1 1
##
## Equal Matrix Constraints:
## ceq Nordstrom Starbucks Microsoft
## Budget -1 -1 -1 -1
## attr(,"na.action")
## Return
## 1
## attr(,"class")
## [1] "omit"
##
## Cardinality Constraints:
## Nordstrom Starbucks Microsoft
## Lower 0 0 0
## Upper 1 1 1
#setSolver(Spec) = "solveRquadprog"
Frontier <- portfolioFrontier(as.timeSeries(ret.ts), Spec, constraints = "Short") #vs "LongOnly" "Short"
## Warning in as.vector(invSigma %*% one)/(one %*% invSigma %*% one): Recycling array of length 1 in vector-array arithmetic is deprecated.
## Use c() or as.vector() instead.
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)
## Cov
## 0.07651345
varRisk(Data, Weights, alpha = 0.05)###VaR of equal weights portfolio
## VaR.5%
## -0.1474633
cvarRisk(Data, Weights, alpha = 0.05)### CVaR for equal weight portfolio
## CVaR.5%
## -0.166304
#===========================================================
# Tangency portfolio with risk free rate = 0.01
#===========================================================
# But the answer seems confusing because it remains the same
# when we change risk free rate to other values!
# ==========================================================
setRiskFreeRate(Spec)<-0.01/12
tgPortfolio <- tangencyPortfolio(ret.ts, Spec, constraints = "Short")
tgPortfolio
##
## Title:
## MV Tangency Portfolio
## Estimator: shrinkEstimator
## Solver: solveRshortExact
## Optimize: minRisk
## Constraints: Short
##
## Portfolio Weights:
## Nordstrom Starbucks Microsoft
## -0.0285 0.2502 0.7783
##
## Covariance Risk Budgets:
## Nordstrom Starbucks Microsoft
## -0.0066 0.1898 0.8168
##
## Target Returns and Risks:
## mean mu Cov Sigma CVaR VaR
## 0.0403 0.0403 0.0879 0.0865 0.1496 0.1146
##
## Description:
## Fri Mar 29 19:11:51 2019 by user:
#=============================================================
# Tangency portfolio by closed form
#==============================================================
Sigma = cov(firm_data1[,2:4])
rf = 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.0418781921
a1 = 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)
##
## Title:
## MV Minimum Variance Portfolio
## Estimator: covEstimator
## Solver: solveRquadprog
## Optimize: minRisk
## Constraints: LongOnly
##
## Portfolio Weights:
## Nordstrom Starbucks Microsoft
## 0.3636 0.1937 0.4427
##
## Covariance Risk Budgets:
## Nordstrom Starbucks Microsoft
## 0.3636 0.1937 0.4427
##
## Target Returns and Risks:
## mean Cov CVaR VaR
## 0.0250 0.0733 0.1288 0.1183
##
## Description:
## Fri Mar 29 19:11:51 2019 by user:
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
}
#==================================================
# Here is the graph output
#===================================================
#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.206 0.279 0.927
lines(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()
#===================================
# no short sales constraint
#===================================
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.ps",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.326 0.297 1.029
lines(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()
#======================================================
# Use systematic investors toolbox (SIT)
# Load Systematic Investor Toolbox (SIT)
# http://www.r-bloggers.com/backtesting-minimum-variance-portfolios/
#=========================================================
#setInternet2(TRUE)
con = gzcon(url('https://github.com/systematicinvestor/SIT/raw/master/sit.gz', 'rb'))
source(con)
close(con)
#*****************************************************************
# Create Constraints
#*****************************************************************
n <- dim(firm_data1.xts)[2]
constraints = new.constraints(n, lb = -Inf, ub = +Inf)
# SUM x.i = 1
constraints = add.constraints(rep(1, n), 1, type = '=', constraints)
ia <- create.historical.ia(firm_data1.xts, 12)
#s0 <- apply(coredata(firm_data1.xts),2,sd)
#ia$cov <- cor(coredata(firm_data1.xts), use='complete.obs',method='pearson') * (s0 %*% t(s0))
weight <- min.risk.portfolio(ia, constraints)
#===================================================
# Try to plot efficient frontier using SIT
# create sample historical input assumptions
ia <- create.historical.ia(firm_data1.xts, 12) # 12 is annual factor for monthly return
# create long-only, fully invested efficient frontier
# 0 <= x.i <= 1
# If short sale allowed: constraints = new.constraints(n, lb = -Inf, ub = +Inf)
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)
# SUM x.i = 1
constraints = add.constraints(rep(1, n), 1, type = '=', constraints)
#
weight <- min.risk.portfolio(ia, constraints)
# create efficient frontier
ifelse(!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':
##
## cross
ef
## $weight
## Nordstrom Starbucks Microsoft
## [1,] 0.36360 0.1937 0.443
## [2,] 0.35417 0.1947 0.451
## [3,] 0.34475 0.1958 0.459
## [4,] 0.33533 0.1969 0.468
## [5,] 0.32590 0.1980 0.476
## [6,] 0.31648 0.1991 0.484
## [7,] 0.30705 0.2002 0.493
## [8,] 0.29763 0.2013 0.501
## [9,] 0.28820 0.2024 0.509
## [10,] 0.27878 0.2035 0.518
## [11,] 0.26935 0.2046 0.526
## [12,] 0.25993 0.2057 0.534
## [13,] 0.25050 0.2068 0.543
## [14,] 0.24108 0.2079 0.551
## [15,] 0.23165 0.2090 0.559
## [16,] 0.22223 0.2101 0.568
## [17,] 0.21280 0.2112 0.576
## [18,] 0.20338 0.2123 0.584
## [19,] 0.19395 0.2134 0.593
## [20,] 0.18453 0.2145 0.601
## [21,] 0.17510 0.2156 0.609
## [22,] 0.16568 0.2167 0.618
## [23,] 0.15625 0.2177 0.626
## [24,] 0.14683 0.2188 0.634
## [25,] 0.13740 0.2199 0.643
## [26,] 0.12798 0.2210 0.651
## [27,] 0.11855 0.2221 0.659
## [28,] 0.10913 0.2232 0.668
## [29,] 0.09970 0.2243 0.676
## [30,] 0.09028 0.2254 0.684
## [31,] 0.08085 0.2265 0.693
## [32,] 0.07143 0.2276 0.701
## [33,] 0.06200 0.2287 0.709
## [34,] 0.05258 0.2298 0.718
## [35,] 0.04315 0.2309 0.726
## [36,] 0.03373 0.2320 0.734
## [37,] 0.02430 0.2331 0.743
## [38,] 0.01488 0.2342 0.751
## [39,] 0.00546 0.2353 0.759
## [40,] 0.00000 0.2264 0.774
## [41,] 0.00000 0.2037 0.796
## [42,] 0.00000 0.1811 0.819
## [43,] 0.00000 0.1585 0.842
## [44,] 0.00000 0.1358 0.864
## [45,] 0.00000 0.1132 0.887
## [46,] 0.00000 0.0905 0.909
## [47,] 0.00000 0.0679 0.932
## [48,] 0.00000 0.0453 0.955
## [49,] 0.00000 0.0226 0.977
## [50,] 0.00000 0.0000 1.000
##
## $return
## [,1]
## [1,] 0.373
## [2,] 0.379
## [3,] 0.384
## [4,] 0.390
## [5,] 0.396
## [6,] 0.401
## [7,] 0.407
## [8,] 0.413
## [9,] 0.418
## [10,] 0.424
## [11,] 0.430
## [12,] 0.436
## [13,] 0.441
## [14,] 0.447
## [15,] 0.453
## [16,] 0.458
## [17,] 0.464
## [18,] 0.470
## [19,] 0.475
## [20,] 0.481
## [21,] 0.487
## [22,] 0.492
## [23,] 0.498
## [24,] 0.504
## [25,] 0.510
## [26,] 0.515
## [27,] 0.521
## [28,] 0.527
## [29,] 0.532
## [30,] 0.538
## [31,] 0.544
## [32,] 0.549
## [33,] 0.555
## [34,] 0.561
## [35,] 0.566
## [36,] 0.572
## [37,] 0.578
## [38,] 0.584
## [39,] 0.589
## [40,] 0.595
## [41,] 0.601
## [42,] 0.606
## [43,] 0.612
## [44,] 0.618
## [45,] 0.623
## [46,] 0.629
## [47,] 0.635
## [48,] 0.640
## [49,] 0.646
## [50,] 0.652
##
## $risk
## [1] 0.254 0.254 0.254 0.254 0.254 0.255 0.255 0.256 0.256 0.257 0.257
## [12] 0.258 0.259 0.259 0.260 0.261 0.262 0.263 0.264 0.266 0.267 0.268
## [23] 0.269 0.271 0.272 0.274 0.275 0.277 0.278 0.280 0.282 0.284 0.286
## [34] 0.287 0.289 0.291 0.293 0.295 0.298 0.300 0.302 0.306 0.309 0.314
## [45] 0.318 0.324 0.329 0.335 0.342 0.349
##
## $name
## [1] "Efficient Frontier"
#====================================================
# David Ruppert's example in his textbook
#===================================================
ifelse(!require(Ecdat), install.packages("Ecdat"), library(Ecdat))
## Loading required package: Ecdat
## Loading required package: Ecfun
##
## Attaching package: 'Ecfun'
## The following object is masked from 'package:base':
##
## sign
##
## Attaching package: 'Ecdat'
## The following object is masked from 'package:datasets':
##
## Orange
## [1] "Ecdat"
ifelse(!require(Ecfun), install.packages("Ecfun"), library(Ecfun))
## [1] "Ecdat"
ifelse(!require(quadprog), install.packages("quadprog"), library(quadprog))
## [1] "Ecdat"
#
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)
## [1] "mts" "ts"
CRSP.df = as.data.frame(CRSPday)
head(CRSP.df)
## year month day ge ibm mobil crsp
## 1 1989 1 3 -0.01676 0.00000 -0.00275 -0.00762
## 2 1989 1 4 0.01705 0.00513 0.00551 0.01302
## 3 1989 1 5 -0.00279 -0.00204 0.00548 0.00281
## 4 1989 1 6 0.00000 -0.00613 0.00272 0.00306
## 5 1989 1 9 0.00000 0.00411 0.00543 0.00163
## 6 1989 1 10 -0.00560 -0.00717 0.00811 -0.00199
#
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
## [1] 0.5512 0.0844 0.3645
lines(c(0,2),mufree+c(0,2)*(muP[ind]-mufree)/sdP[ind],lwd=4,lty=2)