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
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()
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.
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)
##                     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.32826
str(data)
## An 'xts' object on 2007-01-03 08:00:00/2019-03-27 08:00:00 containing:
##   Data: num [1:3079, 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
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.32826
tail(data)
##            JWN.Adjusted SBUX.Adjusted MSFT.Adjusted
## 2019-03-20        43.45         71.63        117.52
## 2019-03-21        43.76         72.26        120.22
## 2019-03-22        42.93         71.96        117.05
## 2019-03-25        43.54         72.30        117.66
## 2019-03-26        43.73         72.96        117.91
## 2019-03-27        44.65         72.74        116.77
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<-
## 
## Attaching package: 'fBasics'
## The following object is masked from 'package:TTR':
## 
##     volatility
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
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
return = firm_data1[,2:4]
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
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)

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]
ret.ts<- timeSeries(return, date1)
chart.CumReturns(ret.ts, legend.loc = 'topleft', main = '')

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:  
##  Portfolio Weights:         NA
##  Target Return:             0.0242391
##  Target Risk:               NA
##  Risk-Free Rate:            0
##  Number of Frontier Points: 50
## 
## Optim List:  
##  Solver:                    solveRshortExact
##  Objective:                 portfolioObjective portfolioReturn portfolioRisk
##  Trace:                     FALSE
setEstimator(Spec)="covMcdEstimator"
setEstimator(Spec)= "covOGKEstimator"
setEstimator(Spec)= "shrinkEstimator"
Spec
## 
## Model List:  
##  Type:                      MV
##  Optimize:                  minRisk
##  Estimator:                 shrinkEstimator
##  Params:                    alpha = 0.05
## 
## Portfolio List:  
##  Portfolio Weights:         NA
##  Target Return:             0.0242391
##  Target Risk:               NA
##  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)
## 
## Title:
##  MV Efficient Portfolio 
##  Estimator:         shrinkEstimator 
##  Solver:            solveRshortExact 
##  Optimize:          minRisk 
##  Constraints:       Short 
## 
## Portfolio Weights:
## Nordstrom Starbucks Microsoft 
##    0.3780    0.2044    0.4176 
## 
## Covariance Risk Budgets:
## Nordstrom Starbucks Microsoft 
##    0.3872    0.2126    0.4001 
## 
## Target Returns and Risks:
##   mean     mu    Cov  Sigma   CVaR    VaR 
## 0.0242 0.0242 0.0734 0.0667 0.1324 0.1243 
## 
## Description:
##  Thu Mar 28 05:03:04 2019 by user: Nomio
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:
##  Thu Mar 28 05:03:05 2019 by user: Nomio
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:
##  Thu Mar 28 05:03:05 2019 by user: Nomio
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
## Return  0.0242391  0.001545085  0.02846068  0.04271153
## Budget -1.0000000 -1.000000000 -1.00000000 -1.00000000
## 
## Cardinality Constraints:
##       Nordstrom Starbucks Microsoft
## Lower         0         0         0
## Upper         1         1         1
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)

weightsPlot(Frontier)

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
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:
##  Thu Mar 28 05:03:09 2019 by user: Nomio
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
ret.tp  = sum(mr.mtx*tp)
ret.tp
## [1] 0.04721981
std.tp = sqrt((t(tp)%*%Sigma)%*%tp)
std.tp
##           [,1]
## [1,] 0.1015897
sharpe.tp = (ret.tp - rf)/std.tp
sharpe.tp
##          [,1]
## [1,] 0.456606
globminSpec <- portfolioSpec()
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:
##  Thu Mar 28 05:03:11 2019 by user: Nomio

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 %”) ```

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
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
}
pdf("3firmExample.pdf",width=6,height=5)  
plot(sdP,muP,type="l",xlim=c(0,0.25),ylim=c(0,0.08),lty=1)
mufree = 0.01/12 # input value of risk-free interest rate
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
ind2 = (sdP == min(sdP)) # find the minimum variance portfolio
ind3 = (muP > muP[ind2])
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
muP1 = seq(min(mu)+.0001,max(mu)-.0001,length=300)
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
}

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)
par(new=TRUE)
## Warning in par(new = TRUE): calling par(new=TRUE) with no plot
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
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)

ind2 = (sdP == min(sdP)) # find the minimum variance portfolio
ind3 = (muP > muP[ind2])
con = gzcon(url('https://github.com/systematicinvestor/SIT/raw/master/sit.gz', 'rb'))
source(con)
close(con)
n <- dim(firm_data1.xts)[2]
constraints = new.constraints(n, lb = -Inf, ub = +Inf)
constraints = add.constraints(rep(1, n), 1, type = '=', constraints)  
ia <- create.historical.ia(firm_data1.xts, 12)
weight <- min.risk.portfolio(ia, constraints)
ia <- create.historical.ia(firm_data1.xts, 12)
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)
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"
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)
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