Nothing
demo/Ch3_ModelingAssetReturns_demo.R
In PCRA: Companion to Portfolio Construction and Risk Analysis
## 1. First, install the devtools package using the drop-down menu:
## RStudio > Tools > Install Packages > devtools
## 2. Next, install the PCRA package FROM THE RStudio CONSOLE using
## devtools::install_github("robustport/PCRA"), and load it with:
library(PCRA)
# NOTE: This PCRA is the latest version. DO NOT INSTALL PCRA from CRAN!
## 3. Install all packages that are arguments of library(packageName) below with:
## RStudio > Tools > Install Packages > PackageName
library(data.table)
library(xts)
library(PerformanceAnalytics)
library(PortfolioAnalytics)
library(foreach)
library(CVXR)
library(RPESE)
library(RPEIF)
library(ggplot2)
library(dplyr)
library(RobStatTM)
## 4. Finally, ensure that you have installed optimalRhoPsi
## FROM THE RStudio CONSOLE using
## devtools::install_github("kjellpk/optimalRhoPsi"), and load it with:
library(optimalRhoPsi)
# NOTE: When loading packages, you can safely ignore the various comments
# that appear in the console in red after each package is loaded, including,
# but not limited to:
# "method overwritten", "using 6 threads", "object is masked", etc. etc.
# YOU ARE NOW READY TO RUN THE CH 2 REPRODUCIBILITY CODE
## Table 2.1
# Largecaps
stockItems <- c("Date", "TickerLast", "CapGroupLast", "Return")
dateRange <- c("1993-01-31", "2015-12-31")
returns <- selectCRSPandSPGMI("monthly", dateRange = dateRange,
stockItems = stockItems, factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "LargeCap",
outputType = "xts")
ret <- coredata(returns) # Package xts
n <- ncol(ret)
acfLarge <- rep(0, n)
for(i in 1:n){
acfLarge[i] <- acf(ret[,i], lag.max = 1, plot = FALSE)$acf[2]
}
muLarge <- mean(acfLarge)
sdLarge <- sd(acfLarge)
# Midcaps
returns <- selectCRSPandSPGMI("monthly",dateRange = dateRange,
stockItems = stockItems, factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "MidCap",
outputType = "xts")
ret <- coredata(returns) # Package xts
n <- ncol(ret)
acfMid <- rep(0, n)
for(i in 1:n){
acfMid[i] <- acf(ret[,i], lag.max = 1,
plot = FALSE)$acf[2]
}
muMid <- mean(acfMid)
sdMid <- sd(acfMid)
## Smallcaps
returns <- selectCRSPandSPGMI("monthly", dateRange = dateRange,
stockItems = stockItems, factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "SmallCap",
outputType = "xts")
ret <- coredata(returns) # Package xts
n <- ncol(ret)
acfSmall <- rep(0, n)
for(i in 1:n){
acfSmall[i] <- acf(ret[,i], lag.max = 1,
plot = FALSE)$acf[2]
}
muSmall <- mean(acfSmall)
sdSmall <- sd(acfSmall)
## Microcaps
returns <- selectCRSPandSPGMI("monthly", dateRange = dateRange,
stockItems = stockItems, factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "MicroCap",
outputType = "xts")
ret <- coredata(returns) # Package xts
n <- ncol(ret)
acfMicro <- rep(0, n)
for(i in 1:n){
acfMicro[i] <- acf(ret[,i], lag.max = 1, plot = FALSE)$acf[2]
}
muMicro <- mean(acfMicro)
sdMicro <- sd(acfMicro)
dat <- cbind(c(muLarge, sdLarge), c(muMid, sdMid),
c(muSmall, sdSmall), c(muMicro, sdMicro))
dat <- round(dat,3)
dat <- data.frame(dat)
names(dat) <- c("LargeCap", "MidCap", "SmallCap", "MicroCap")
row.names(dat) <- c(" Mean Lag-1 Acf", "StdDev Lag-1 Acf")
dat
## Figure 2.1
dat <- list(acfLarge, acfMid, acfSmall, acfMicro)
names(dat) <- c("87 LargeCaps", "67 MidCaps", "106 SmallCaps", "34 MicroCaps")
boxplot(dat, varwidth = TRUE, col = "cyan")
## Figure 2.2
returns <- PerformanceAnalytics::edhec
returns <- returns["2002-01-31/2019-12-31", -13]
names(returns) <- c("CA", "CTA", "DIST", "EM", "EMN", "ED",
"FIA", "GM", "LSE", "MA", "RV", "SS")
PCRA::tsPlotMP(returns)
# range(index(returns))
Ret <- coredata(returns)
n <- ncol(Ret)
acfRet <- rep(0, n)
for(i in 1:n){
acfRet[i] <- acf(Ret[,i], lag.max = 1,
plot = FALSE)$acf[2]
}
hist(acfRet, main = "EDHEC Hedge Fund Indexes",
xlab = "Lag-1 ACF Values")
## Figure 2.3
# names(acfRet) <- names(returns)
# (names(sort(acfRet)[1:4]))
# delete "CTA", "GM", "SS", "EMN"
returns8 <- returns[ , c("CA", "DIST", "EM", "ED",
"FIA", "LSE", "MA", "RV")]
PCRA::tsPlotMP(returns8, yname = "RETURNS",
stripText.cex = 0.7, axis.cex = 0.7)
## Figure 2.4
acf(returns8$EM, main = "EM", lag.max = 5)
## Figure 2.5
muVol <- c(.20, .10, .15, .04)
wts <- seq(0,1, .01)
efront2Asset <- function(wts, rho, muVol = c(.20, .10, .15, .04))
{
sigma1 <- muVol[1]
mu1 <- muVol[2]
sigma2 <- muVol[3]
mu2 <- muVol[4]
n <- length(wts)
efront <- data.frame(matrix(rep(0, 3*n), ncol = 3))
names(efront) <- c("SIGMA", "MU", "WTS")
w <- wts
for(i in 1:n){
mu <- w[i]*mu1 + (1 - w[i])*mu2
var <- w[i]^2*sigma1^2 + 2*w[i]*(1 - w[i])*rho*sigma1*sigma2 + (1 - w[i])^2*sigma2^2
sigma <- sqrt(var)
efront[i,] <- c(sigma, mu, w[i])
}
return(efront)
}
ef <- efront2Asset(wts, 0, muVol = muVol)
gmv <- ef[ef$SIGMA == min(ef$SIGMA),]
xlab <- expression(sigma [P])
ylab <- expression(mu [P])
par(pty = "s")
plot(ef$SIGMA, ef$MU, type = "l", xlab = xlab, ylab = ylab,
xlim = c(0, .25), ylim = c(0.03, .11), lwd = 2, cex.lab = 1.5)
points(muVol[c(1, 3)], muVol[c(2, 4)], pch = 19, cex = 1.3)
points(gmv, pch = 19, cex = 1.3)
text(.04, .10, expression(paste(rho, " = 0")), cex = 1.5)
text(0.12, .0616, adj = c(1, NA), "MinRisk ", cex = 1.1)
text(0.13, .0616, adj = c(0,NA), "(.12, .0616)",cex = 1.1)
text(0.2, .1, adj = c(0, NA), " (.20, .10)", cex = 1.1)
text(0.15, .04, adj = c(0, NA)," (.15, .04)", cex = 1.1)
## Figure 2.6
muVol <- c(.20, .10, .15, .04)
wts <- seq(0, 1, .01)
ef <- efront2Asset(wts, 0, muVol = muVol)
ef1 <- efront2Asset(wts, 1, muVol = muVol)
ef2 <- efront2Asset(wts, -1, muVol = muVol)
gmv <- ef[ef$SIGMA == min(ef$SIGMA),]
gmv2 <- ef2[ef2$SIGMA == min(ef2$SIGMA),]
xlab <- expression(sigma [P])
ylab <- expression(mu [P])
par(pty = "s")
plot(ef$SIGMA, ef$MU, type = "l", xlab = xlab, ylab = ylab,
xlim = c(0, .25), ylim = c(0.03, .11), lwd = 2, cex.lab = 1.5)
points(muVol[c(1, 3)], muVol[c(2, 4)], pch = 19, cex = 1.5)
points(gmv, pch = 19, cex = 1.3)
text(0.2, .1, adj = c(0, NA)," (.20, .10)",cex = 1.1)
text(0.15, .04, adj = c(0,NA)," (.15, .04)", cex = 1.1)
lines(ef1$SIGMA, ef1$MU, lty = 2, lwd = 2)
lines(ef2$SIGMA, ef2$MU, lty = 2, lwd = 2)
points(gmv2, pch = 19, cex = 1.3)
text(.12, .07, expression(paste(rho, " = 0 ")), adj = c(1, NA), cex = 1.5)
text(.02, .08, expression(paste(rho, " = -1 ")), adj = c(0, NA), cex = 1.5)
text(.18, .07, expression(paste(rho, " = +1 ")), adj = c(0, NA), cex = 1.5)
## Figure 2.7
muVol <- c(.20, .10, .15, .04)
wts <- seq(0, 1, .01)
efLO <- efront2Asset(wts, 0, muVol = muVol)
wts <- seq(1, 1.25, .01)
efSS <- efront2Asset(wts, 0,muVol = muVol)
gmv <- ef[ef$SIGMA == min(ef$SIGMA),]
maxMu <- ef[ef$MU == max(ef$MU),]
maxMuSS <- efSS[efSS$MU == max(efSS$MU),]
xlab <- expression(sigma [P])
ylab <- expression(mu [P])
par(pty = "s")
plot(efLO$SIGMA, efLO$MU, type = "l", xlab = xlab, ylab = ylab,
xlim=c(0, .40), ylim=c(.02, .13), lwd = 2, cex.lab = 1.5)
lines(efSS$SIGMA, efSS$MU, lty = "dashed", lwd = 2)
points(gmv[1:2], pch = 19, cex = 1.3)
points(maxMu[1:2], pch = 19, cex = 1.3)
points(maxMuSS[1:2], pch = 19, cex = 1.3)
text(.04, .12, expression(paste(rho, " = 0")), cex = 1.5)
text(gmv[1:2], adj = c(0, NA),
paste(" (",toString(round(gmv[1:2],2)),")"), cex = 1.1)
text(maxMu[1:2], adj = c(0, NA),
paste(" (",toString(maxMu[1:2]),")"), cex = 1.1)
text(maxMuSS[1:2], adj = c(0, NA),
paste(" (",toString(round(maxMuSS[1:2], 2)),")"), cex = 1.1)
## Figure 2.8
volMu1 <- c(.20, .10)
volMu2 <- c(.15, .04)
volMu3 <- c(.10, .02)
names(volMu1) <- c("SIGMA", "MU")
names(volMu2) <- c("SIGMA", "MU")
names(volMu3) <- c("SIGMA", "MU")
wts <- seq(0, 1, .01)
ef1 <- efront2Asset(wts, 0, muVol = c(volMu1, volMu2))
ef2 <- efront2Asset(wts, 0, muVol = c(volMu1, volMu3))
ef3 <- efront2Asset(wts, 0, muVol = c(volMu2, volMu3))
xlab <- expression(sigma [P])
ylab <- expression(mu [P])
par(pty = "s")
plot(ef1$SIGMA,ef1$MU,type = "l", xlab = xlab, ylab = ylab,
xlim=c(0,.25), ylim=c(0,.11), lwd = 2, cex.lab = 1.5)
lines(ef2$SIGMA, ef2$MU, lty = 2, lwd = 2.0)
lines(ef3$SIGMA, ef3$MU, lty = 3, lwd = 2.0)
xy <- rbind(volMu1, volMu2, volMu3)
points(xy, pch = 19, cex = 1.3)
text(volMu1[1] + 0.01, volMu1[2], adj = c(0,NA), toString(volMu1), cex = 1.1)
text(volMu2[1] + 0.01, volMu2[2], adj = c(0,NA), toString(volMu2), cex = 1.1)
text(volMu3[1] + 0.01, volMu3[2], adj = c(0,NA), toString(volMu3), cex = 1.1)
text(.04, .10, expression(paste(rho, " = 0")), cex = 1.5)
## Figure 2.9
volMu1 <- c(.20,.10)
volMu2 <- c(.15,.04)
volMu3 <- c(.10,.02)
vol <- c(volMu1[1], volMu2[1], volMu3[1])
mu <- c(volMu1[2], volMu2[2], volMu3[2])
corrMat0 <- matrix(rep(0, 9),nrow = 3) + diag(rep(1, 3))
covMat0 <- diag(vol) %*% corrMat0 %*% diag(vol)
n <- 500
port <- matrix(rep(0, 2*n), ncol = 2)
dimnames(port)[[2]] = c("SIG.P", "MU.P")
wts = hitandrun::simplex.sample(3, n)$samples
for(i in 1:n) {
x <- wts[i,]
port[i, 1] <- sqrt(x%*%covMat0%*%x)
port[i, 1]
port[i, 2] <- x%*%mu
}
xlab <- expression(sigma [P])
ylab <- expression(mu [P])
plot(port[, 1], port[, 2], xlim = c(0, .25), ylim = c(0, .11),
xlab = xlab, ylab = ylab, pch = 20, cex = .7, cex.lab = 1.5)
points(vol,mu, pch = 19, cex = 1.3)
text(volMu1[1] + 0.01, volMu1[2], adj = c(0,NA), toString(volMu1), cex = 1.1)
text(volMu2[1] + 0.01, volMu2[2], adj = c(0,NA), toString(volMu2), cex = 1.1)
text(volMu3[1] + 0.01, volMu3[2], adj = c(0,NA), toString(volMu3), cex = 1.1)
text(.04, .10, expression(paste(rho, " = 0")), cex = 1.5)
## Example 2.3
muRet <- c(.10, .04, .02)
volRet <- c(.20, .15, .10)
corrRet <- diag(c(1, 1, 1))
PCRA::mathGmvMuCov(muRet,volRet,corrRet,digits = 3)
## Figure 2.10
# There is currently no code for this Figure, and while
# IBM and XOM are in our CRSP data, GE is not.
# Will consider replacing GE with a stock in the CRSP data.
## Example 2.6 Figures 2.8 - 2.10 with the following steps
# Get xts object of 106 smallcap stocks, and the Market ("MktIndexCRSP") in
# stocksCRSP for 1997 - 2010, and use the third group of 10 of these
# to compute Gmv portfolios. Change name "MktIndexCRSP" to "Market".
## Figure 2.11
library(data.table)
stockItems <- c("Date","TickerLast","CapGroupLast","Return","MktIndexCRSP")
dateRange <- c("1997-01-31","2010-12-31")
stocksDat <- PCRA::selectCRSPandSPGMI("monthly",dateRange = dateRange,
stockItems = stockItems,
factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "SmallCap",
outputType = "xts")
returns10Mkt <- stocksDat[, c(21:30,107)]
names(returns10Mkt)[11] <- "Market"
tsPlotMP(returns10Mkt,scaleType = "free",layout = c(2, 6),stripText.cex = .45,
axis.cex = 0.4, lwd = 0.5)
## Figures 2.12 and 2.13
# Use PortfolioAnalytics, and related functions, to compute and plot
# time series of: (1) GmvLS portfolio weights wtsGmvLS, (2) Combined
# GmvLS and Market returns, (3) cumulative gross returns of GmvLS and
# Market portfolios.
# Create GmvLS portfolio specs
returns <- returns10Mkt[, 1:10]
Market <- returns10Mkt[, 11]
funds <- colnames(returns)
# The following portfolio specification functions and optimization
# function are from the PortfolioAnalytics package
pspec <- portfolio.spec(assets=funds)
pspec.fi <- add.constraint(pspec, type="full_investment")
pspec.gmvLS <- add.objective(pspec.fi, type="risk", name="var")
# Optimize Portfolio at Monthly Rebalancing and 5-Year Training
bt.gmvLS <- optimize.portfolio.rebalancing(returns, pspec.gmvLS,
optimize_method = "CVXR",
rebalance_on = "months",
training_period = 60,
rolling_window = 60,
trace = TRUE)
# Extract time series of portfolio weights
wtsGmvLS <- extractWeights(bt.gmvLS)
# Compute rebalancing GmvLS arithmetic returns
# The Return.rebalancing function is from PerformanceAnalytics
GmvLS <- Return.rebalancing(returns, wtsGmvLS)
# Combine GmvLS and Market returns and plot their time series
ret.comb <- na.omit(merge.xts(GmvLS, Market, all=F))
names(ret.comb) <- c("GmvLS", "Market")
# Figure 2.12
tsPlotMP(wtsGmvLS, layout = c(2,5), scaleType = "same",
stripText.cex = 0.7, axis.cex = .7)
# Figure 2.13
tsPlotMP(ret.comb, scaleType = "same", stripText.cex = .7, axis.cex = .7)
## Figure 2.14
# Compute cumulative gross portfolio returns
R <- ret.comb
geometric <- TRUE
c.xts <- if ( geometric ) {
cumprod(1+R)
} else {
1 + cumsum(R)
}
# Plot cumulative gross returns of GmvLS and Market portfolios
# Original code contributed by Peter Carl
p <- plot.xts(c.xts[,1], col="black", main = "Cumulative Returns",
grid.ticks.lwd=1, grid.ticks.lty = "dotted", grid.ticks.on = "years",
labels.col="grey20", cex.axis=0.8, format.labels = "%b\n%Y",
ylim = c(min(c.xts), max(c.xts)))
p <- addSeries(c.xts[,2], on=1, lwd=2, col="darkred", lty="dashed")
p <- addLegend("topleft", on = 1,
legend.names = names(c.xts),
lty = c(1,2), lwd = rep(2, NCOL(c.xts)),
col = c("black", "darkred"),
bty = "o", box.col = "white",
bg=rgb(t(col2rgb("white")), alpha = 200,
maxColorValue = 255) )
d.xts <- PerformanceAnalytics::Drawdowns(R)
p <- xts::addSeries(d.xts[,1], col="darkblue", lwd=2, main="Drawdown",
ylim = c(min(d.xts), 0) )
p <- xts::addSeries(d.xts[,2], on=2, lwd=2, col="darkred", lty="dashed")
# panel 1 and 2 ylim
## ylim1 <- c(p$Env$ylim[[2]][1], p$Env$ylim[[2]][2]) No longer works
## ylim2 <- c(p$Env$ylim[[4]][1], p$Env$ylim[[4]][2]) No longer works
ylim1 <- p$Env$panels[[1]]$ylim
ylim2 <- p$Env$panels[[2]]$ylim
ylim <- c(ylim1, ylim2)
# get longest drawdown dates for xts object
dt <- table.Drawdowns(R, top = 1) # just want to find the worst drawdown
dt2 <- t(dt[, c("From", "To")])
x <- as.vector(dt2[, NCOL(dt2)])
y <- as.xts(matrix(rep(ylim, length(x)), ncol=length(ylim), byrow=TRUE),
order.by=as.Date(x))
i=1
p <- xts::addPolygon(y[i:(i+1), 1:2], on=-1, col="lightgrey") # top panel
p <- xts::addPolygon(y[i:(i+1), 3:4], on=-2, col="lightgrey") # lower panel
p
## Table 2.2
dt2mat <- table.Drawdowns(R, top = 2) # find the worst two drawdowns
dt2mat[,4] <- round(dt2mat[,4], 2)
dt2 <- data.frame(dt2mat)[, 1:5]
names(dt2) <- c("Begin", "Minimum", "End", "Depth", "Months")
dt2 <- dt2[, c(1:3, 5, 4)]
dt2
## Figure 2.15
levg <- levgLongShort(wtsGmvLS)
plot.zoo(levg, ylim = c(0.0,1.5), ylab = "Leverage")
abline(h = 1.0, lty = "dotted")
## Figure 2.16 Left-Hand Plot
GmvLS.TO <- 100*turnOver(wtsGmvLS)
plot.zoo(GmvLS.TO, ylim = c(0,60), ylab = "TURNOVER (%)",
xlab = "", cex.axis = 1.5, cex.lab = 1.5)
abline(h = mean(GmvLS.TO), lty = "dashed")
text(as.Date("2004-01-31"), 50, "Mean Turnover = 13.4 (%)", cex = 1.5)
## Figure 2.16 Right-Hand Plot
GmvLS.DIV <- 100*divHHI(wtsGmvLS)
plot.zoo(GmvLS.DIV,ylim = c(0,100),lwd = 1.5, ylab = "DIV(%)",
xlab = "", cex.axis = 1.5, cex.lab = 1.5)
abline(h = mean(GmvLS.DIV),lty = "dashed")
text(as.Date("2006-01-31"), 90, "Mean Diversification = 63.5 (%)", cex = 1.5)
## Table 2.3
# Risk & Performance Estimator Standard Errors package
SD12 <- SD.SE(ret.comb, se.method = "IFiid")
SD12 <- printSE(SD12, round.digit = 4)
SSD12 <- SemiSD.SE(ret.comb, se.method = "IFiid")
SSD12 <- printSE(SSD12, round.digit = 4)
ES12 <- ES.SE(ret.comb, se.method = "IFiid")
ES12 <- printSE(ES12, round.digit = 4)
VaR12 <- VaR.SE(ret.comb, se.method = "IFiid")
VaR12 <- printSE(VaR12, round.digit = 4)
# VaR12[,1] <- -VaR12[,1]
RM <- 100*rbind(SD12,SSD12,ES12,VaR12)
colnames(RM) <- c("Estimate (%)","StdError (%)")
rownames(RM) <- c("GmvLS SD", "Market SD",
"GmvLS SSD", "Market SSD",
"GmvLS ES", "Market ES",
"GmvLS VaR", "Market VaR")
RM <- as.data.frame(RM)
RM
## Figure 2.17
# Risk-free rates were not negligible before 2009
stockItems <- c("Date", "TickerLast", "Return", "Ret13WkBill")
returnsAll <- selectCRSPandSPGMI("monthly", stockItems = stockItems,
factorItems = NULL, outputType = "xts")
riskFree <- returnsAll[ , "Ret13WkBill"]
tsPlotMP(riskFree, yname = "RISK-FREE RATE")
## Figure 2.18
# Variation of risk-free rate for 1995 through 2000
returns <- returnsAll["1995-01-31/2000-12-31"]
x <- sort(apply(returns, 2, mean))
x0 <- x[x <= 0.007 & x >= 0.005] # Results in 21 stocks & choose FMC
ret <- returns[, c("FMC", "Ret13WkBill")]
names(ret)[2] <- "Risk-Free"
tsPlotMP(ret, yname = "RETURNS", scaleType = "free")
## Table 2.4
SR12 <- SR.SE(ret.comb, se.method = "IFiid")
SR12 <- printSE(SR12, round.digit = 2)
DSR12 <- DSR.SE(ret.comb, se.method = "IFiid")
DSR12 <- printSE(DSR12, round.digit = 2)
SoR12 <- SoR.SE(ret.comb, se.method = "IFiid")
SoR12 <- printSE(SoR12, round.digit = 2)
ESratio12 <- ESratio.SE(ret.comb, se.method = "IFiid")
ESratio12 <- printSE(ESratio12, round.digit = 2)
Ratios <- rbind(SR12,DSR12,SoR12,ESratio12)
colnames(Ratios) <- c("Estimate","StdError")
rownames(Ratios) <- c("GmvLS SR", "Market SR",
"GmvLS DSR","Market DSR",
"GmvLS SOR","Market SOR",
"GmvLS ESR","Market ESR")
Ratios <- as.data.frame(Ratios)
Ratios
## Figure 2.19
data(edhec, package = "PerformanceAnalytics")
colnames(edhec) <- c("CA", "CTAG", "DIS", "EM", "EMN", "ED", "FIA",
"GM", "LS", "MA", "RV", "SS", "FoF")
par(mfrow = c(1, 2))
outSD <- IF.SD(returns = edhec$CA, evalShape = T, IFplot = T)
outSR <- IF.SemiSD(returns = edhec$CA, evalShape = T, IFplot = T)
par(mfrow = c(1, 1))
## Figure 2.20
muRet <- c(.10, .04, .02)
volRet <- c(.20, .15, .10)
corrRet <- diag(c(1, 1, 1))
mathEfrontRiskyMuCov(muRet, volRet, corrRet, efront.only = F)
text(0.07, 0.095, "EFFICIENT FRONTIER", cex = 1.2)
arrows(0.07, 0.09, .10, .06, length = 0.1, lwd= 1.0)
## Table 2.5
muRet <- c(.10, .04, .02)
volRet <- c(.20, .15, .10)
corrRet <- diag(c(1, 1, 1))
efront <- mathEfrontRiskyMuCov(muRet, volRet, corrRet, npoints = 5,
values = T, display = F)
mu.efront <- efront$mu.efront
wtsEfront <- mathWtsEfrontRiskyMuCov(muRet, volRet, corrRet,
mu.efront, digits = 3)
wtsEfront
## Figure 2.21
nColor <- 4
barplotWts(as.matrix(wtsEfront), legend.text = T, ylab = "WEIGHTS",
col = topo.colors(nColor), bar.ylim = c(-1, 2),
cex.lab = 1.2, cex.axis = 1.3)
## Figure 2.22
returns10 <- returns10Mkt[,-11]
efront <- mathEfrontRisky(returns10, display = T, cexGmv = 1.0,
cexPoints = 1.1, cexText = 0.9)
## Figure 2.23
efront10 <- mathEfrontRisky(returns10,npoints = 5, display = F, values = TRUE)
mu.efront <- efront10$mu.efront
wtsEfront <- mathWtsEfrontRisky(returns10, mu.efront,digits = 3)
barplotWts(as.matrix(wtsEfront), legend.text = T, ylab = "WEIGHTS",
col = topo.colors(10), bar.ylim = c(-1.5,3.0),cex.lab = 1.2,
cex.axis = 1.3)
## Figure 2.24
stockItems <- c("Date", "TickerLast", "CapGroupLast", "Return", "Ret13WkBill")
dateRange <- c("1997-01-31", "2010-12-31")
stocksDat <- selectCRSPandSPGMI("monthly", dateRange = dateRange,
stockItems = stockItems,
factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "SmallCap",
outputType = "xts")
returns10andRF <- stocksDat[, c(21:30,107)]
names(returns10andRF)[11] <- "RiskFree"
tsPlotMP(returns10andRF,scaleType = "free",layout = c(2,6),stripText.cex = .45,
axis.cex = 0.4,lwd = 0.5)
## Figure 2.25
rf <- mean(returns10andRF[, 11])
returns10 <- returns10andRF[, -11]
mathEfrontCashRisky(returns10, rf = rf, cexPoints = 1.0)
## Figure 2.26
rf <- mean(returns10andRF[, 11])
returns10 <- returns10andRF[, -11]
wtsEfront = mathEfrontCashRisky(returns10, plot.efront = FALSE, values = TRUE)
barplotWts(as.matrix(wtsEfront),legend.text = T,col = topo.colors(3),
ylab = "Weights",xlab = "VOL", bar.ylim = c(-0.5, 1.5))
## Calculation of risk aversion and risk tolerance values for Figure 2.25,
## reported at end of paragraph below (2.154)
rf <- .005
C <- var(returns10)
mu.stocks <- apply(returns10, 2, mean)
mue <- mu.stocks - rf
a <- solve(C, mue)
lambda <- sum(a)
lambda # Risk aversion value
1/lambda # Risk tolerance value
## Figure 2.27
plot(FRBinterestRates, xaxt = "n", xlab = "", ylab = "InterestRates (%)")
axis(side = 1,at = seq(1930,2020,by=10), labels = seq(1930, 2020, by=10))
grid()
## Figure 2.28
rf <- mean(returns10andRF[,11])
returns10 <- returns10andRF[,-11]
mathEfront(returns10, rf = rf, mu.max = .035, sigma.max = .19, cexText = 0.8, npoints = 100)
## Figure 2.29
rf <- 0.03
rf_lend <- 0.04
rf_borrow <- 0.06
er_port <- 0.07
leverage <- seq(0, 2, .1)
er_rf <- rf + leverage * (er_port - rf)
er_1 <- rf_lend + leverage * (er_port - rf_lend)
er_2 <- rf_borrow + leverage * (er_port - rf_borrow)
df <- data.frame("Leverage" = leverage,
"Single Risk Free Rate for Borrowing and Lending" = er_rf,
"Different Risk Free Rates for Borrowing and Lending" = pmin(er_1, er_2))
df_melt <- reshape2::melt(df, id.vars = "Leverage", variable.name = "Risk_Free_Rate")
df_melt[["Risk_Free_Rate"]] <- gsub("\\.", " ", df_melt[["Risk_Free_Rate"]])
ggplot(df_melt, aes(x = Leverage, y = value )) +
geom_line(aes(color = Risk_Free_Rate, linetype = Risk_Free_Rate), linewidth = 1) +
labs(x = "Leverage", y = "Expected Return",
color = "Risk Free Rate", linetype = "Risk Free Rate") +
theme_bw() +
theme(legend.position = "inside",
legend.position.inside = c(0.4, 0.9),
legend.title = element_text(size = 16),
legend.text = element_text(size = 14),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14))
## Figure 2.30
volMu1 <- c(.20, .10)
volMu2 <- c(.15, .04)
volMu3 <- c(.10, .02)
vol <- c(volMu1[1], volMu2[1], volMu3[1])
mu <- c(volMu1[2], volMu2[2], volMu3[2])
corrMat <- matrix(rep(0, 9), nrow = 3)+ diag(rep(1, 3))
V <- diag(vol) %*% corrMat %*% diag(vol)
# Compute IR for the three stocks
one <- rep(1, nrow(V))
z1 <- solve(V, one) # Vinv*one
z2 <- solve(V, mu) # Vinv*mu
a <- as.numeric(t(mu) %*% z1) # a = mu*Vinv*one
b <- as.numeric(t(mu) %*% z2) # b = mu*Vinv*mu
cc <- as.numeric(t(one) %*% z1) # c = one*Vinv*one
d <- b * cc - a^2
muGmv <- a/cc
sigmaGmv <- 1/sqrt(cc)
IR <- sigmaGmv*sqrt(d)
# Plot active efficient frontier
sigmaA <- seq(0, 20, 1)
muA <- IR*sigmaA
xlab <- "Active Volatility (%)"
ylab <- "Active Mean Return (%)"
plot(sigmaA, muA, xlim = c(0, 21), ylim = c(0, 8.5), type = "l", lwd = 1.5,
xlab = xlab, ylab = ylab, xaxs = "i", yaxs = "i", cex.lab = 1.3)
text(2, 7, pos = 4, "IR = slope of line = .36", cex = 1.5)
## Table 2.6
wGmv <- z1/cc
w1 <- z2/a
mu1 <- b/a
sigmaA <- c(.02,.05,.10)
wA <- as.matrix((IR*sigmaA/(mu1 - muGmv))%*%t(w1-wGmv))
rowSum <- apply(wA,1,sum)
wA <- cbind(wA,rowSum)
wA <- round(wA,3)
wA.df <- data.frame(wtA1=wA[,1],wtA2=wA[,2],wtA3=wA[,3],
wtAsum=wA[,4])
rnames <- c("TE 2% ", "TE 5% ", "TE 10% ")
row.names(wA.df) <- rnames
wA.df
## Figure 2.31
# Actively managed frontier dominated by efficient frontier
volB <- 0.12
muB <- 0.045
varGmv <- sigmaGmv^2
muGmv <- muGmv
mu1 <- mu1
varB <- volB^2
muA <- seq(-.02, 0.06, .001)
sigmaA <- muA/IR
muPA <- muB + muA
varA <- sigmaA^2
const <- (2/(mu1 - muGmv)) * varGmv * (muB/muGmv - 1)
varPA <- varB + varA + muA*const
sigmaPA <- sqrt(varPA)
# Plot using mathEfrontRiskyMuCov for the efficient frontier
mathEfrontRiskyMuCov(mu, vol, corrMat, efront.only = T, display = T)
lines(sigmaPA, muPA, type = "l", lwd = 1.5)
points(volB, muB, pch = 16, cex = 1.3)
text(volB, muB, "B", pos = 4, cex = 1.3)
text(0.163, 0.06, "ACTIVELY MANAGED", pos = 4, cex = 1.1)
arrows(0.165, 0.06, sigmaPA[50], muPA[50], length = .1, lwd = 1.5)
## Figure 2.32
# The two plots in this Figure are from Jorion(2003)
# Proper citation will be added.
## Figure 2.33
# These two plots were created by the PCRA authors
## Figure 2.34
# Left-Hand Plot Code
powerUtilityPlots <- function()
{
x <- seq(.01,3,.01)
y <- log(x)
lwd <- 1.0
plot(x, y, axes=F, type = "l", ylim =c(-8, 2), lwd = lwd, xlab = "v", ylab = "U(v)")
axis(side = 1, pos = 0)
axis(side = 2, pos = 0)
gamma <- -0.5
shift <- 1
y <- (x^gamma - shift)/gamma
lines(x, y, lty = 8, lwd = lwd)
gamma <- 0.5
y = (x^gamma - shift)/gamma
lines(x, y, lty = 3, lwd = lwd)
abline(v = 0)
legend(1.2, -5.5, c("Gamma = .5", "Log Utility","Gamma = -.5"),
lty = c(3, 1, 8), lwd = 1.0)
}
powerUtilityPlots()
# Right-Hand Plot Code
quadraticUtilityPlot <- function()
{
v <- seq(0, 1.5, .01)
u <- v - v^2
ylim <- c(-0.7, 0.4)
plot(v, u, type = "l", ylim = ylim, xlab = "v", ylab = "U(v)", lwd = 1.5)
abline(v = .50, lty = "dotted")
abline(h = .25, lty = "dotted")
}
quadraticUtilityPlot()
## Figure 2.35
# Left-Hand Plot Code
rm1 <- 0.18
Beta1 <- c(1.53, 1.36, 1.24, 1.17, 1.06, 0.92, 0.84, 0.76, 0.63, 0.48)
Mu1 <- c(0.26, 0.22, 0.21, 0.21, 0.18, 0.17, 0.16, 0.15, 0.13, 0.12)
Sigma1 <- c(0.49, 0.43, 0.39, 0.37, 0.33, 0.29, 0.27, 0.24, 0.20, 0.17)
SML1 <- rm1 * Beta1
df1 <- data.frame( Beta1, Mu1, Sigma1, SML1)
df1a <- data.frame(x=1, y = rm1)
p1 <- ggplot(df1) +
geom_point(aes(x = Beta1, y = Mu1), color = "black") +
geom_line(aes( x = Beta1, y = SML1), color = "gray20") +
labs(x = "Beta", y = "Mean Excess Return", title = "1931 \u2013 1965") +
theme(plot.title = element_text(hjust = 0.5)) +
annotate(geom="text", x=1.19, y=0.15,
label="CAPM Security Market Line", color="gray20") +
annotate(geom="text", x=.78, y=0.185,
label="Market Portfolio", color="dodgerblue4")
p1 + geom_point(data = df1a, aes(x = x, y = y),
color = "dodgerblue4", shape = 17, size = 3)
# Rigt-Hand Plot Code
rm2 <- 0.08
Beta2 <- c(1.50, 1.30, 1.17, 1.09, 1.03, 0.95, 0.87, 0.78, 0.67, 0.51)
Mu2 <- c(0.06, 0.08, 0.08, 0.08, 0.08, 0.08, 0.07, 0.07, 0.07, 0.06)
Sigma2 <- c(0.31, 0.26, 0.24, 0.22, 0.21, 0.19, 0.18, 0.16, 0.14, 0.12)
SML2 <- rm2 * Beta2
df2 <- data.frame(Beta2, Mu2, Sigma2, SML2)
df2a <- data.frame(x = 1, y = rm2)
p2 <- ggplot(df2) +
geom_point(aes(x = Beta2, y = Mu2), color = "black") +
geom_line(aes( x = Beta2, y = SML2), color = "gray20") +
ylim(0, 0.12) +
labs(x = "Beta", y = "Mean Excess Return", title = "1965 \u2013 1991") +
theme(plot.title = element_text(hjust = 0.5)) +
annotate(geom = "text", x = 1.12, y = 0.055,
label = "CAPM Security Market Line", color = "gray20") +
annotate(geom = "text", x = .80, y = 0.095,
label = "Market Portfolio", color = "dodgerblue4")
p2 + geom_point(data = df2a, aes(x = x, y = y),
color = "dodgerblue4", shape = 17, size = 3)
## Figure 2.36
# Left-Hand Plot Code
ggplot(df1, aes(x = Beta1, y = Sigma1)) +
geom_point() +
geom_smooth(formula = 'y ~ x', method='lm', se = FALSE, linewidth = 0.6, color = "gray20") +
labs(x = "Beta", y = "Standard Deviation") +
geom_text(x = 1.05, y = 0.27, label = "sigma %~~% 0.32 ~ beta", parse=TRUE)
# Right-Hand Plot Code
ggplot(df2, aes(x = Beta2, y = Sigma2)) +
geom_point() +
geom_smooth(formula = 'y ~ x', method='lm', se = FALSE, linewidth = 0.6, color = "gray20") +
labs(x = "Beta", y = "Standard Deviation") +
geom_text(x = 1.05, y = 0.18, label = "sigma %~~% 0.2 ~ beta", parse=TRUE)
## Table 2.7
df <- data.frame(matrix(" ", nrow = 6, ncol = 4))
df$X1 <- c("1/31--12/39", "1/40--12/49", "1/50--12/59", "1/60--12/69", "1/70--12/79","1/80--12/91")
df$X2 <- c(-0.05, 0.03, 0.08, 0.03, 0.01, 0.09)
df$X3 <- c(0.17, 0.10, 0.06, 0.07, 0.10, 0.08)
df$X4 <- c(-0.94, 1.06, 4.25, 1.32, 0.18, 3.90)
#Rename rows and columns and reformat the table
# colnames(df) <- c("Period", "$\\mu_{e}$", "$\\sigma_{e}$", "$t( \\mu )$") # In the text
colnames(df) <- c("Period", "Mean(Excess Return)", "Std. Dev(Excess Return)", "t(Mean(Excess Return))")
df
## Figure 2.37
stocksCRSPweekly <- getPCRAData("stocksCRSPweekly")
dateRange <- c("2004-01-01", "2005-12-31")
stockItems <- c("Date", "TickerLast", "CapGroupLast",
"Return", "MktIndexCRSP", "Ret13WkBill")
returnsAll <- selectCRSPandSPGMI("weekly",
dateRange = dateRange,
stockItems = stockItems,
factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "SmallCap",
outputType = "xts")
returns <- returnsAll[ , 1:10]
tsPlotMP(returns, scaleType = "free",layout = c(2,5),
stripText.cex = .45, axis.cex = 0.4,lwd = 0.5)
## Figure 2.38
pspec <- portfolio.spec(assets = names(returns))
pspecFI <- add.constraint(pspec, type = "full_investment")
pspecLO <- add.constraint(portfolio = pspecFI, type = "long_only")
pspecESLO5pct <- add.objective(pspecLO, type = "risk", name = "ES",
arguments = list(p = 0.050))
# Increase n.portfolios below for more accurate vertical dot-dash line
lty <- c("dashed", "solid", "dotted", "dotdash")
col <- c("red", "black", "darkgreen", "darkgreen")
chart.EfficientFrontierCompare(returns, pspecESLO5pct, risk_type = "ES",
guideline = TRUE, cex.axis = 1.2,
match.col = c("StdDev", "ES"),
n.portfolios = 10,
lwd = c(1.3, 1.4, 1.3, 1.0),
col = col, lty = lty,
xlim = c(0.02, 0.08), ylim = c(0.0, 0.012),
legend.loc = "topleft", main = NULL)
## Figure 2.39
chart.EfficientFrontierCompare(returns, pspecESLO5pct, risk_type = "StdDev",
guideline = TRUE, cex.axis = 1.2,
match.col = c("ES", "StdDev"),
n.portfolios = 10,
lwd = c(1.3, 1.4, 1.3, 1.0),
col = col, lty = lty,
xlim = c(0.01, 0.06), ylim = c(0.0, 0.012),
legend.loc = "topleft", main = NULL)
## Figure 2.40
x <- seq(-4.9, 4.9, by = 0.001)
ccopt <- computeTuningPsi_modOpt(0.95)
plot(x, wgt_modOpt(x, cc = ccopt), type = "l",
xlab = "x", ylab = "", cex = 1.5, cex.lab = 1.5)
## Figure 2.41
data(edhec)
hfnames <- c("CA", "CTA", "DIS", "EM", "EMN", "ED", "FIA",
"GM", "LSE", "MA", "RV", "SS", "FOF")
names(edhec) <- hfnames
retLongFIA <- edhec[, "FIA"]
retFIA <- retLongFIA['1998-01-31/1999-12-31', ]
index(retFIA) <- as.yearmon(index(retFIA))
# Plot FIA returns with sample mean and robust mean
mu <- 100*mean(retFIA)
se.mu <- 100*sd(retFIA)/sqrt(24)
x <- RobStatTM::locScaleM(retFIA, eff = .95)
muRob <- 100*x$mu
se.muRob <- 100*x$std.mu
plot.zoo(retFIA, type ="b", xlab = "", ylab = "FIA Returns")
abline(h = muRob/100, col = "blue")
abline(h = mu/100, lty = "dashed", col ="red")
legend(1999.2, -.03, legend = c("Robust Mean", "Sample Mean"), lwd = c(1, 2),
lty = c("solid", "dashed"), col = c("blue", "red"), bty = "n", cex = 1.3)
## Table 2.8
tstat.mu <- mu/se.mu
tstat.muRob <- muRob/se.muRob
SR.classic <- tstat.mu/sqrt(24)
SR.Rob <- tstat.muRob/sqrt(24)
row1 <- round(c(mu, se.mu, tstat.mu, SR.classic), 2)
row2 <- round(c(muRob, se.muRob, tstat.muRob, SR.Rob), 2)
meanEsts <- data.frame(rbind(row1, row2))
names(meanEsts) <- c("Estimate (%)", "Std. Error (%)", "t-Stat", "Sharpe Ratio")
row.names(meanEsts) <- c("Sample Mean", "Robust Mean")
meanEsts
## Figure 2.42
# mOpt 95% Efficiency Tuning Constant
ccModOpt <- computeTuningPsi_modOpt(0.95)
# mOpt M-scale weight function
wgt_mOptScale <- function(x) {
rho_modOpt(x, cc = ccModOpt)/(x^2)}
wgt_mOptScaleMax <- wgt_mOptScale(0.0001)
# Plot mOptScale Weight Function
x <- seq(-5.5, 5.5, 0.01)
ylim <- c(0, 1.4)
ylab <- "w_scale,mOpt(x)"
plot(x,wgt_mOptScale(x)/wgt_mOptScaleMax, ylim = ylim,
ylab = ylab, type = "l", cex.lab = 1.1)
abline(h = 1.0, lty = "dotted")
## Figure 2.43
data(edhec, package = "PerformanceAnalytics")
hfnames <- c("CA", "CTA", "DIS", "EM", "EMN", "ED", "FIA",
"GM", "LSE", "MA", "RV", "SS", "FOF")
names(edhec) <- hfnames
range(index(edhec))
edhec <- edhec[ , 1:12]
returns <- 100*edhec['1998-01-31/1999-12-31']
# Plot edhec returns for 1998-1999
tsPlotMP(returns, type = "l", stripText.cex = 0.7, axis.cex = 0.7)
## Table 2.9
StdDev <- apply(returns,2,sd)
MADM <- apply(returns,2,mad)
resid <- returns - median(returns)
RobSD <- apply(resid, 2, scaleM, family = "mopt")
SDestsMat <- cbind(StdDev, MADM, RobSD)
SDestsMat <- round(SDestsMat,2)
SDests <- data.frame(hfnames[1:12], SDestsMat)
names(SDests) <- c("HFindex", "StdDev", "MADM", "RobSD")
row.names(SDests) <- NULL
SDests
Try the PCRA package in your browser
Any scripts or data that you put into this service are public.
PCRA documentation built on July 15, 2026, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## 1. First, install the devtools package using the drop-down menu:
## RStudio > Tools > Install Packages > devtools
## 2. Next, install the PCRA package FROM THE RStudio CONSOLE using
## devtools::install_github("robustport/PCRA"), and load it with:
library(PCRA)
# NOTE: This PCRA is the latest version. DO NOT INSTALL PCRA from CRAN!
## 3. Install all packages that are arguments of library(packageName) below with:
## RStudio > Tools > Install Packages > PackageName
library(data.table)
library(xts)
library(PerformanceAnalytics)
library(PortfolioAnalytics)
library(foreach)
library(CVXR)
library(RPESE)
library(RPEIF)
library(ggplot2)
library(dplyr)
library(RobStatTM)
## 4. Finally, ensure that you have installed optimalRhoPsi
## FROM THE RStudio CONSOLE using
## devtools::install_github("kjellpk/optimalRhoPsi"), and load it with:
library(optimalRhoPsi)
# NOTE: When loading packages, you can safely ignore the various comments
# that appear in the console in red after each package is loaded, including,
# but not limited to:
# "method overwritten", "using 6 threads", "object is masked", etc. etc.
# YOU ARE NOW READY TO RUN THE CH 2 REPRODUCIBILITY CODE
## Table 2.1
# Largecaps
stockItems <- c("Date", "TickerLast", "CapGroupLast", "Return")
dateRange <- c("1993-01-31", "2015-12-31")
returns <- selectCRSPandSPGMI("monthly", dateRange = dateRange,
stockItems = stockItems, factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "LargeCap",
outputType = "xts")
ret <- coredata(returns) # Package xts
n <- ncol(ret)
acfLarge <- rep(0, n)
for(i in 1:n){
acfLarge[i] <- acf(ret[,i], lag.max = 1, plot = FALSE)$acf[2]
}
muLarge <- mean(acfLarge)
sdLarge <- sd(acfLarge)
# Midcaps
returns <- selectCRSPandSPGMI("monthly",dateRange = dateRange,
stockItems = stockItems, factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "MidCap",
outputType = "xts")
ret <- coredata(returns) # Package xts
n <- ncol(ret)
acfMid <- rep(0, n)
for(i in 1:n){
acfMid[i] <- acf(ret[,i], lag.max = 1,
plot = FALSE)$acf[2]
}
muMid <- mean(acfMid)
sdMid <- sd(acfMid)
## Smallcaps
returns <- selectCRSPandSPGMI("monthly", dateRange = dateRange,
stockItems = stockItems, factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "SmallCap",
outputType = "xts")
ret <- coredata(returns) # Package xts
n <- ncol(ret)
acfSmall <- rep(0, n)
for(i in 1:n){
acfSmall[i] <- acf(ret[,i], lag.max = 1,
plot = FALSE)$acf[2]
}
muSmall <- mean(acfSmall)
sdSmall <- sd(acfSmall)
## Microcaps
returns <- selectCRSPandSPGMI("monthly", dateRange = dateRange,
stockItems = stockItems, factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "MicroCap",
outputType = "xts")
ret <- coredata(returns) # Package xts
n <- ncol(ret)
acfMicro <- rep(0, n)
for(i in 1:n){
acfMicro[i] <- acf(ret[,i], lag.max = 1, plot = FALSE)$acf[2]
}
muMicro <- mean(acfMicro)
sdMicro <- sd(acfMicro)
dat <- cbind(c(muLarge, sdLarge), c(muMid, sdMid),
c(muSmall, sdSmall), c(muMicro, sdMicro))
dat <- round(dat,3)
dat <- data.frame(dat)
names(dat) <- c("LargeCap", "MidCap", "SmallCap", "MicroCap")
row.names(dat) <- c(" Mean Lag-1 Acf", "StdDev Lag-1 Acf")
dat
## Figure 2.1
dat <- list(acfLarge, acfMid, acfSmall, acfMicro)
names(dat) <- c("87 LargeCaps", "67 MidCaps", "106 SmallCaps", "34 MicroCaps")
boxplot(dat, varwidth = TRUE, col = "cyan")
## Figure 2.2
returns <- PerformanceAnalytics::edhec
returns <- returns["2002-01-31/2019-12-31", -13]
names(returns) <- c("CA", "CTA", "DIST", "EM", "EMN", "ED",
"FIA", "GM", "LSE", "MA", "RV", "SS")
PCRA::tsPlotMP(returns)
# range(index(returns))
Ret <- coredata(returns)
n <- ncol(Ret)
acfRet <- rep(0, n)
for(i in 1:n){
acfRet[i] <- acf(Ret[,i], lag.max = 1,
plot = FALSE)$acf[2]
}
hist(acfRet, main = "EDHEC Hedge Fund Indexes",
xlab = "Lag-1 ACF Values")
## Figure 2.3
# names(acfRet) <- names(returns)
# (names(sort(acfRet)[1:4]))
# delete "CTA", "GM", "SS", "EMN"
returns8 <- returns[ , c("CA", "DIST", "EM", "ED",
"FIA", "LSE", "MA", "RV")]
PCRA::tsPlotMP(returns8, yname = "RETURNS",
stripText.cex = 0.7, axis.cex = 0.7)
## Figure 2.4
acf(returns8$EM, main = "EM", lag.max = 5)
## Figure 2.5
muVol <- c(.20, .10, .15, .04)
wts <- seq(0,1, .01)
efront2Asset <- function(wts, rho, muVol = c(.20, .10, .15, .04))
{
sigma1 <- muVol[1]
mu1 <- muVol[2]
sigma2 <- muVol[3]
mu2 <- muVol[4]
n <- length(wts)
efront <- data.frame(matrix(rep(0, 3*n), ncol = 3))
names(efront) <- c("SIGMA", "MU", "WTS")
w <- wts
for(i in 1:n){
mu <- w[i]*mu1 + (1 - w[i])*mu2
var <- w[i]^2*sigma1^2 + 2*w[i]*(1 - w[i])*rho*sigma1*sigma2 + (1 - w[i])^2*sigma2^2
sigma <- sqrt(var)
efront[i,] <- c(sigma, mu, w[i])
}
return(efront)
}
ef <- efront2Asset(wts, 0, muVol = muVol)
gmv <- ef[ef$SIGMA == min(ef$SIGMA),]
xlab <- expression(sigma [P])
ylab <- expression(mu [P])
par(pty = "s")
plot(ef$SIGMA, ef$MU, type = "l", xlab = xlab, ylab = ylab,
xlim = c(0, .25), ylim = c(0.03, .11), lwd = 2, cex.lab = 1.5)
points(muVol[c(1, 3)], muVol[c(2, 4)], pch = 19, cex = 1.3)
points(gmv, pch = 19, cex = 1.3)
text(.04, .10, expression(paste(rho, " = 0")), cex = 1.5)
text(0.12, .0616, adj = c(1, NA), "MinRisk ", cex = 1.1)
text(0.13, .0616, adj = c(0,NA), "(.12, .0616)",cex = 1.1)
text(0.2, .1, adj = c(0, NA), " (.20, .10)", cex = 1.1)
text(0.15, .04, adj = c(0, NA)," (.15, .04)", cex = 1.1)
## Figure 2.6
muVol <- c(.20, .10, .15, .04)
wts <- seq(0, 1, .01)
ef <- efront2Asset(wts, 0, muVol = muVol)
ef1 <- efront2Asset(wts, 1, muVol = muVol)
ef2 <- efront2Asset(wts, -1, muVol = muVol)
gmv <- ef[ef$SIGMA == min(ef$SIGMA),]
gmv2 <- ef2[ef2$SIGMA == min(ef2$SIGMA),]
xlab <- expression(sigma [P])
ylab <- expression(mu [P])
par(pty = "s")
plot(ef$SIGMA, ef$MU, type = "l", xlab = xlab, ylab = ylab,
xlim = c(0, .25), ylim = c(0.03, .11), lwd = 2, cex.lab = 1.5)
points(muVol[c(1, 3)], muVol[c(2, 4)], pch = 19, cex = 1.5)
points(gmv, pch = 19, cex = 1.3)
text(0.2, .1, adj = c(0, NA)," (.20, .10)",cex = 1.1)
text(0.15, .04, adj = c(0,NA)," (.15, .04)", cex = 1.1)
lines(ef1$SIGMA, ef1$MU, lty = 2, lwd = 2)
lines(ef2$SIGMA, ef2$MU, lty = 2, lwd = 2)
points(gmv2, pch = 19, cex = 1.3)
text(.12, .07, expression(paste(rho, " = 0 ")), adj = c(1, NA), cex = 1.5)
text(.02, .08, expression(paste(rho, " = -1 ")), adj = c(0, NA), cex = 1.5)
text(.18, .07, expression(paste(rho, " = +1 ")), adj = c(0, NA), cex = 1.5)
## Figure 2.7
muVol <- c(.20, .10, .15, .04)
wts <- seq(0, 1, .01)
efLO <- efront2Asset(wts, 0, muVol = muVol)
wts <- seq(1, 1.25, .01)
efSS <- efront2Asset(wts, 0,muVol = muVol)
gmv <- ef[ef$SIGMA == min(ef$SIGMA),]
maxMu <- ef[ef$MU == max(ef$MU),]
maxMuSS <- efSS[efSS$MU == max(efSS$MU),]
xlab <- expression(sigma [P])
ylab <- expression(mu [P])
par(pty = "s")
plot(efLO$SIGMA, efLO$MU, type = "l", xlab = xlab, ylab = ylab,
xlim=c(0, .40), ylim=c(.02, .13), lwd = 2, cex.lab = 1.5)
lines(efSS$SIGMA, efSS$MU, lty = "dashed", lwd = 2)
points(gmv[1:2], pch = 19, cex = 1.3)
points(maxMu[1:2], pch = 19, cex = 1.3)
points(maxMuSS[1:2], pch = 19, cex = 1.3)
text(.04, .12, expression(paste(rho, " = 0")), cex = 1.5)
text(gmv[1:2], adj = c(0, NA),
paste(" (",toString(round(gmv[1:2],2)),")"), cex = 1.1)
text(maxMu[1:2], adj = c(0, NA),
paste(" (",toString(maxMu[1:2]),")"), cex = 1.1)
text(maxMuSS[1:2], adj = c(0, NA),
paste(" (",toString(round(maxMuSS[1:2], 2)),")"), cex = 1.1)
## Figure 2.8
volMu1 <- c(.20, .10)
volMu2 <- c(.15, .04)
volMu3 <- c(.10, .02)
names(volMu1) <- c("SIGMA", "MU")
names(volMu2) <- c("SIGMA", "MU")
names(volMu3) <- c("SIGMA", "MU")
wts <- seq(0, 1, .01)
ef1 <- efront2Asset(wts, 0, muVol = c(volMu1, volMu2))
ef2 <- efront2Asset(wts, 0, muVol = c(volMu1, volMu3))
ef3 <- efront2Asset(wts, 0, muVol = c(volMu2, volMu3))
xlab <- expression(sigma [P])
ylab <- expression(mu [P])
par(pty = "s")
plot(ef1$SIGMA,ef1$MU,type = "l", xlab = xlab, ylab = ylab,
xlim=c(0,.25), ylim=c(0,.11), lwd = 2, cex.lab = 1.5)
lines(ef2$SIGMA, ef2$MU, lty = 2, lwd = 2.0)
lines(ef3$SIGMA, ef3$MU, lty = 3, lwd = 2.0)
xy <- rbind(volMu1, volMu2, volMu3)
points(xy, pch = 19, cex = 1.3)
text(volMu1[1] + 0.01, volMu1[2], adj = c(0,NA), toString(volMu1), cex = 1.1)
text(volMu2[1] + 0.01, volMu2[2], adj = c(0,NA), toString(volMu2), cex = 1.1)
text(volMu3[1] + 0.01, volMu3[2], adj = c(0,NA), toString(volMu3), cex = 1.1)
text(.04, .10, expression(paste(rho, " = 0")), cex = 1.5)
## Figure 2.9
volMu1 <- c(.20,.10)
volMu2 <- c(.15,.04)
volMu3 <- c(.10,.02)
vol <- c(volMu1[1], volMu2[1], volMu3[1])
mu <- c(volMu1[2], volMu2[2], volMu3[2])
corrMat0 <- matrix(rep(0, 9),nrow = 3) + diag(rep(1, 3))
covMat0 <- diag(vol) %*% corrMat0 %*% diag(vol)
n <- 500
port <- matrix(rep(0, 2*n), ncol = 2)
dimnames(port)[[2]] = c("SIG.P", "MU.P")
wts = hitandrun::simplex.sample(3, n)$samples
for(i in 1:n) {
x <- wts[i,]
port[i, 1] <- sqrt(x%*%covMat0%*%x)
port[i, 1]
port[i, 2] <- x%*%mu
}
xlab <- expression(sigma [P])
ylab <- expression(mu [P])
plot(port[, 1], port[, 2], xlim = c(0, .25), ylim = c(0, .11),
xlab = xlab, ylab = ylab, pch = 20, cex = .7, cex.lab = 1.5)
points(vol,mu, pch = 19, cex = 1.3)
text(volMu1[1] + 0.01, volMu1[2], adj = c(0,NA), toString(volMu1), cex = 1.1)
text(volMu2[1] + 0.01, volMu2[2], adj = c(0,NA), toString(volMu2), cex = 1.1)
text(volMu3[1] + 0.01, volMu3[2], adj = c(0,NA), toString(volMu3), cex = 1.1)
text(.04, .10, expression(paste(rho, " = 0")), cex = 1.5)
## Example 2.3
muRet <- c(.10, .04, .02)
volRet <- c(.20, .15, .10)
corrRet <- diag(c(1, 1, 1))
PCRA::mathGmvMuCov(muRet,volRet,corrRet,digits = 3)
## Figure 2.10
# There is currently no code for this Figure, and while
# IBM and XOM are in our CRSP data, GE is not.
# Will consider replacing GE with a stock in the CRSP data.
## Example 2.6 Figures 2.8 - 2.10 with the following steps
# Get xts object of 106 smallcap stocks, and the Market ("MktIndexCRSP") in
# stocksCRSP for 1997 - 2010, and use the third group of 10 of these
# to compute Gmv portfolios. Change name "MktIndexCRSP" to "Market".
## Figure 2.11
library(data.table)
stockItems <- c("Date","TickerLast","CapGroupLast","Return","MktIndexCRSP")
dateRange <- c("1997-01-31","2010-12-31")
stocksDat <- PCRA::selectCRSPandSPGMI("monthly",dateRange = dateRange,
stockItems = stockItems,
factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "SmallCap",
outputType = "xts")
returns10Mkt <- stocksDat[, c(21:30,107)]
names(returns10Mkt)[11] <- "Market"
tsPlotMP(returns10Mkt,scaleType = "free",layout = c(2, 6),stripText.cex = .45,
axis.cex = 0.4, lwd = 0.5)
## Figures 2.12 and 2.13
# Use PortfolioAnalytics, and related functions, to compute and plot
# time series of: (1) GmvLS portfolio weights wtsGmvLS, (2) Combined
# GmvLS and Market returns, (3) cumulative gross returns of GmvLS and
# Market portfolios.
# Create GmvLS portfolio specs
returns <- returns10Mkt[, 1:10]
Market <- returns10Mkt[, 11]
funds <- colnames(returns)
# The following portfolio specification functions and optimization
# function are from the PortfolioAnalytics package
pspec <- portfolio.spec(assets=funds)
pspec.fi <- add.constraint(pspec, type="full_investment")
pspec.gmvLS <- add.objective(pspec.fi, type="risk", name="var")
# Optimize Portfolio at Monthly Rebalancing and 5-Year Training
bt.gmvLS <- optimize.portfolio.rebalancing(returns, pspec.gmvLS,
optimize_method = "CVXR",
rebalance_on = "months",
training_period = 60,
rolling_window = 60,
trace = TRUE)
# Extract time series of portfolio weights
wtsGmvLS <- extractWeights(bt.gmvLS)
# Compute rebalancing GmvLS arithmetic returns
# The Return.rebalancing function is from PerformanceAnalytics
GmvLS <- Return.rebalancing(returns, wtsGmvLS)
# Combine GmvLS and Market returns and plot their time series
ret.comb <- na.omit(merge.xts(GmvLS, Market, all=F))
names(ret.comb) <- c("GmvLS", "Market")
# Figure 2.12
tsPlotMP(wtsGmvLS, layout = c(2,5), scaleType = "same",
stripText.cex = 0.7, axis.cex = .7)
# Figure 2.13
tsPlotMP(ret.comb, scaleType = "same", stripText.cex = .7, axis.cex = .7)
## Figure 2.14
# Compute cumulative gross portfolio returns
R <- ret.comb
geometric <- TRUE
c.xts <- if ( geometric ) {
cumprod(1+R)
} else {
1 + cumsum(R)
}
# Plot cumulative gross returns of GmvLS and Market portfolios
# Original code contributed by Peter Carl
p <- plot.xts(c.xts[,1], col="black", main = "Cumulative Returns",
grid.ticks.lwd=1, grid.ticks.lty = "dotted", grid.ticks.on = "years",
labels.col="grey20", cex.axis=0.8, format.labels = "%b\n%Y",
ylim = c(min(c.xts), max(c.xts)))
p <- addSeries(c.xts[,2], on=1, lwd=2, col="darkred", lty="dashed")
p <- addLegend("topleft", on = 1,
legend.names = names(c.xts),
lty = c(1,2), lwd = rep(2, NCOL(c.xts)),
col = c("black", "darkred"),
bty = "o", box.col = "white",
bg=rgb(t(col2rgb("white")), alpha = 200,
maxColorValue = 255) )
d.xts <- PerformanceAnalytics::Drawdowns(R)
p <- xts::addSeries(d.xts[,1], col="darkblue", lwd=2, main="Drawdown",
ylim = c(min(d.xts), 0) )
p <- xts::addSeries(d.xts[,2], on=2, lwd=2, col="darkred", lty="dashed")
# panel 1 and 2 ylim
## ylim1 <- c(p$Env$ylim[[2]][1], p$Env$ylim[[2]][2]) No longer works
## ylim2 <- c(p$Env$ylim[[4]][1], p$Env$ylim[[4]][2]) No longer works
ylim1 <- p$Env$panels[[1]]$ylim
ylim2 <- p$Env$panels[[2]]$ylim
ylim <- c(ylim1, ylim2)
# get longest drawdown dates for xts object
dt <- table.Drawdowns(R, top = 1) # just want to find the worst drawdown
dt2 <- t(dt[, c("From", "To")])
x <- as.vector(dt2[, NCOL(dt2)])
y <- as.xts(matrix(rep(ylim, length(x)), ncol=length(ylim), byrow=TRUE),
order.by=as.Date(x))
i=1
p <- xts::addPolygon(y[i:(i+1), 1:2], on=-1, col="lightgrey") # top panel
p <- xts::addPolygon(y[i:(i+1), 3:4], on=-2, col="lightgrey") # lower panel
p
## Table 2.2
dt2mat <- table.Drawdowns(R, top = 2) # find the worst two drawdowns
dt2mat[,4] <- round(dt2mat[,4], 2)
dt2 <- data.frame(dt2mat)[, 1:5]
names(dt2) <- c("Begin", "Minimum", "End", "Depth", "Months")
dt2 <- dt2[, c(1:3, 5, 4)]
dt2
## Figure 2.15
levg <- levgLongShort(wtsGmvLS)
plot.zoo(levg, ylim = c(0.0,1.5), ylab = "Leverage")
abline(h = 1.0, lty = "dotted")
## Figure 2.16 Left-Hand Plot
GmvLS.TO <- 100*turnOver(wtsGmvLS)
plot.zoo(GmvLS.TO, ylim = c(0,60), ylab = "TURNOVER (%)",
xlab = "", cex.axis = 1.5, cex.lab = 1.5)
abline(h = mean(GmvLS.TO), lty = "dashed")
text(as.Date("2004-01-31"), 50, "Mean Turnover = 13.4 (%)", cex = 1.5)
## Figure 2.16 Right-Hand Plot
GmvLS.DIV <- 100*divHHI(wtsGmvLS)
plot.zoo(GmvLS.DIV,ylim = c(0,100),lwd = 1.5, ylab = "DIV(%)",
xlab = "", cex.axis = 1.5, cex.lab = 1.5)
abline(h = mean(GmvLS.DIV),lty = "dashed")
text(as.Date("2006-01-31"), 90, "Mean Diversification = 63.5 (%)", cex = 1.5)
## Table 2.3
# Risk & Performance Estimator Standard Errors package
SD12 <- SD.SE(ret.comb, se.method = "IFiid")
SD12 <- printSE(SD12, round.digit = 4)
SSD12 <- SemiSD.SE(ret.comb, se.method = "IFiid")
SSD12 <- printSE(SSD12, round.digit = 4)
ES12 <- ES.SE(ret.comb, se.method = "IFiid")
ES12 <- printSE(ES12, round.digit = 4)
VaR12 <- VaR.SE(ret.comb, se.method = "IFiid")
VaR12 <- printSE(VaR12, round.digit = 4)
# VaR12[,1] <- -VaR12[,1]
RM <- 100*rbind(SD12,SSD12,ES12,VaR12)
colnames(RM) <- c("Estimate (%)","StdError (%)")
rownames(RM) <- c("GmvLS SD", "Market SD",
"GmvLS SSD", "Market SSD",
"GmvLS ES", "Market ES",
"GmvLS VaR", "Market VaR")
RM <- as.data.frame(RM)
RM
## Figure 2.17
# Risk-free rates were not negligible before 2009
stockItems <- c("Date", "TickerLast", "Return", "Ret13WkBill")
returnsAll <- selectCRSPandSPGMI("monthly", stockItems = stockItems,
factorItems = NULL, outputType = "xts")
riskFree <- returnsAll[ , "Ret13WkBill"]
tsPlotMP(riskFree, yname = "RISK-FREE RATE")
## Figure 2.18
# Variation of risk-free rate for 1995 through 2000
returns <- returnsAll["1995-01-31/2000-12-31"]
x <- sort(apply(returns, 2, mean))
x0 <- x[x <= 0.007 & x >= 0.005] # Results in 21 stocks & choose FMC
ret <- returns[, c("FMC", "Ret13WkBill")]
names(ret)[2] <- "Risk-Free"
tsPlotMP(ret, yname = "RETURNS", scaleType = "free")
## Table 2.4
SR12 <- SR.SE(ret.comb, se.method = "IFiid")
SR12 <- printSE(SR12, round.digit = 2)
DSR12 <- DSR.SE(ret.comb, se.method = "IFiid")
DSR12 <- printSE(DSR12, round.digit = 2)
SoR12 <- SoR.SE(ret.comb, se.method = "IFiid")
SoR12 <- printSE(SoR12, round.digit = 2)
ESratio12 <- ESratio.SE(ret.comb, se.method = "IFiid")
ESratio12 <- printSE(ESratio12, round.digit = 2)
Ratios <- rbind(SR12,DSR12,SoR12,ESratio12)
colnames(Ratios) <- c("Estimate","StdError")
rownames(Ratios) <- c("GmvLS SR", "Market SR",
"GmvLS DSR","Market DSR",
"GmvLS SOR","Market SOR",
"GmvLS ESR","Market ESR")
Ratios <- as.data.frame(Ratios)
Ratios
## Figure 2.19
data(edhec, package = "PerformanceAnalytics")
colnames(edhec) <- c("CA", "CTAG", "DIS", "EM", "EMN", "ED", "FIA",
"GM", "LS", "MA", "RV", "SS", "FoF")
par(mfrow = c(1, 2))
outSD <- IF.SD(returns = edhec$CA, evalShape = T, IFplot = T)
outSR <- IF.SemiSD(returns = edhec$CA, evalShape = T, IFplot = T)
par(mfrow = c(1, 1))
## Figure 2.20
muRet <- c(.10, .04, .02)
volRet <- c(.20, .15, .10)
corrRet <- diag(c(1, 1, 1))
mathEfrontRiskyMuCov(muRet, volRet, corrRet, efront.only = F)
text(0.07, 0.095, "EFFICIENT FRONTIER", cex = 1.2)
arrows(0.07, 0.09, .10, .06, length = 0.1, lwd= 1.0)
## Table 2.5
muRet <- c(.10, .04, .02)
volRet <- c(.20, .15, .10)
corrRet <- diag(c(1, 1, 1))
efront <- mathEfrontRiskyMuCov(muRet, volRet, corrRet, npoints = 5,
values = T, display = F)
mu.efront <- efront$mu.efront
wtsEfront <- mathWtsEfrontRiskyMuCov(muRet, volRet, corrRet,
mu.efront, digits = 3)
wtsEfront
## Figure 2.21
nColor <- 4
barplotWts(as.matrix(wtsEfront), legend.text = T, ylab = "WEIGHTS",
col = topo.colors(nColor), bar.ylim = c(-1, 2),
cex.lab = 1.2, cex.axis = 1.3)
## Figure 2.22
returns10 <- returns10Mkt[,-11]
efront <- mathEfrontRisky(returns10, display = T, cexGmv = 1.0,
cexPoints = 1.1, cexText = 0.9)
## Figure 2.23
efront10 <- mathEfrontRisky(returns10,npoints = 5, display = F, values = TRUE)
mu.efront <- efront10$mu.efront
wtsEfront <- mathWtsEfrontRisky(returns10, mu.efront,digits = 3)
barplotWts(as.matrix(wtsEfront), legend.text = T, ylab = "WEIGHTS",
col = topo.colors(10), bar.ylim = c(-1.5,3.0),cex.lab = 1.2,
cex.axis = 1.3)
## Figure 2.24
stockItems <- c("Date", "TickerLast", "CapGroupLast", "Return", "Ret13WkBill")
dateRange <- c("1997-01-31", "2010-12-31")
stocksDat <- selectCRSPandSPGMI("monthly", dateRange = dateRange,
stockItems = stockItems,
factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "SmallCap",
outputType = "xts")
returns10andRF <- stocksDat[, c(21:30,107)]
names(returns10andRF)[11] <- "RiskFree"
tsPlotMP(returns10andRF,scaleType = "free",layout = c(2,6),stripText.cex = .45,
axis.cex = 0.4,lwd = 0.5)
## Figure 2.25
rf <- mean(returns10andRF[, 11])
returns10 <- returns10andRF[, -11]
mathEfrontCashRisky(returns10, rf = rf, cexPoints = 1.0)
## Figure 2.26
rf <- mean(returns10andRF[, 11])
returns10 <- returns10andRF[, -11]
wtsEfront = mathEfrontCashRisky(returns10, plot.efront = FALSE, values = TRUE)
barplotWts(as.matrix(wtsEfront),legend.text = T,col = topo.colors(3),
ylab = "Weights",xlab = "VOL", bar.ylim = c(-0.5, 1.5))
## Calculation of risk aversion and risk tolerance values for Figure 2.25,
## reported at end of paragraph below (2.154)
rf <- .005
C <- var(returns10)
mu.stocks <- apply(returns10, 2, mean)
mue <- mu.stocks - rf
a <- solve(C, mue)
lambda <- sum(a)
lambda # Risk aversion value
1/lambda # Risk tolerance value
## Figure 2.27
plot(FRBinterestRates, xaxt = "n", xlab = "", ylab = "InterestRates (%)")
axis(side = 1,at = seq(1930,2020,by=10), labels = seq(1930, 2020, by=10))
grid()
## Figure 2.28
rf <- mean(returns10andRF[,11])
returns10 <- returns10andRF[,-11]
mathEfront(returns10, rf = rf, mu.max = .035, sigma.max = .19, cexText = 0.8, npoints = 100)
## Figure 2.29
rf <- 0.03
rf_lend <- 0.04
rf_borrow <- 0.06
er_port <- 0.07
leverage <- seq(0, 2, .1)
er_rf <- rf + leverage * (er_port - rf)
er_1 <- rf_lend + leverage * (er_port - rf_lend)
er_2 <- rf_borrow + leverage * (er_port - rf_borrow)
df <- data.frame("Leverage" = leverage,
"Single Risk Free Rate for Borrowing and Lending" = er_rf,
"Different Risk Free Rates for Borrowing and Lending" = pmin(er_1, er_2))
df_melt <- reshape2::melt(df, id.vars = "Leverage", variable.name = "Risk_Free_Rate")
df_melt[["Risk_Free_Rate"]] <- gsub("\\.", " ", df_melt[["Risk_Free_Rate"]])
ggplot(df_melt, aes(x = Leverage, y = value )) +
geom_line(aes(color = Risk_Free_Rate, linetype = Risk_Free_Rate), linewidth = 1) +
labs(x = "Leverage", y = "Expected Return",
color = "Risk Free Rate", linetype = "Risk Free Rate") +
theme_bw() +
theme(legend.position = "inside",
legend.position.inside = c(0.4, 0.9),
legend.title = element_text(size = 16),
legend.text = element_text(size = 14),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14))
## Figure 2.30
volMu1 <- c(.20, .10)
volMu2 <- c(.15, .04)
volMu3 <- c(.10, .02)
vol <- c(volMu1[1], volMu2[1], volMu3[1])
mu <- c(volMu1[2], volMu2[2], volMu3[2])
corrMat <- matrix(rep(0, 9), nrow = 3)+ diag(rep(1, 3))
V <- diag(vol) %*% corrMat %*% diag(vol)
# Compute IR for the three stocks
one <- rep(1, nrow(V))
z1 <- solve(V, one) # Vinv*one
z2 <- solve(V, mu) # Vinv*mu
a <- as.numeric(t(mu) %*% z1) # a = mu*Vinv*one
b <- as.numeric(t(mu) %*% z2) # b = mu*Vinv*mu
cc <- as.numeric(t(one) %*% z1) # c = one*Vinv*one
d <- b * cc - a^2
muGmv <- a/cc
sigmaGmv <- 1/sqrt(cc)
IR <- sigmaGmv*sqrt(d)
# Plot active efficient frontier
sigmaA <- seq(0, 20, 1)
muA <- IR*sigmaA
xlab <- "Active Volatility (%)"
ylab <- "Active Mean Return (%)"
plot(sigmaA, muA, xlim = c(0, 21), ylim = c(0, 8.5), type = "l", lwd = 1.5,
xlab = xlab, ylab = ylab, xaxs = "i", yaxs = "i", cex.lab = 1.3)
text(2, 7, pos = 4, "IR = slope of line = .36", cex = 1.5)
## Table 2.6
wGmv <- z1/cc
w1 <- z2/a
mu1 <- b/a
sigmaA <- c(.02,.05,.10)
wA <- as.matrix((IR*sigmaA/(mu1 - muGmv))%*%t(w1-wGmv))
rowSum <- apply(wA,1,sum)
wA <- cbind(wA,rowSum)
wA <- round(wA,3)
wA.df <- data.frame(wtA1=wA[,1],wtA2=wA[,2],wtA3=wA[,3],
wtAsum=wA[,4])
rnames <- c("TE 2% ", "TE 5% ", "TE 10% ")
row.names(wA.df) <- rnames
wA.df
## Figure 2.31
# Actively managed frontier dominated by efficient frontier
volB <- 0.12
muB <- 0.045
varGmv <- sigmaGmv^2
muGmv <- muGmv
mu1 <- mu1
varB <- volB^2
muA <- seq(-.02, 0.06, .001)
sigmaA <- muA/IR
muPA <- muB + muA
varA <- sigmaA^2
const <- (2/(mu1 - muGmv)) * varGmv * (muB/muGmv - 1)
varPA <- varB + varA + muA*const
sigmaPA <- sqrt(varPA)
# Plot using mathEfrontRiskyMuCov for the efficient frontier
mathEfrontRiskyMuCov(mu, vol, corrMat, efront.only = T, display = T)
lines(sigmaPA, muPA, type = "l", lwd = 1.5)
points(volB, muB, pch = 16, cex = 1.3)
text(volB, muB, "B", pos = 4, cex = 1.3)
text(0.163, 0.06, "ACTIVELY MANAGED", pos = 4, cex = 1.1)
arrows(0.165, 0.06, sigmaPA[50], muPA[50], length = .1, lwd = 1.5)
## Figure 2.32
# The two plots in this Figure are from Jorion(2003)
# Proper citation will be added.
## Figure 2.33
# These two plots were created by the PCRA authors
## Figure 2.34
# Left-Hand Plot Code
powerUtilityPlots <- function()
{
x <- seq(.01,3,.01)
y <- log(x)
lwd <- 1.0
plot(x, y, axes=F, type = "l", ylim =c(-8, 2), lwd = lwd, xlab = "v", ylab = "U(v)")
axis(side = 1, pos = 0)
axis(side = 2, pos = 0)
gamma <- -0.5
shift <- 1
y <- (x^gamma - shift)/gamma
lines(x, y, lty = 8, lwd = lwd)
gamma <- 0.5
y = (x^gamma - shift)/gamma
lines(x, y, lty = 3, lwd = lwd)
abline(v = 0)
legend(1.2, -5.5, c("Gamma = .5", "Log Utility","Gamma = -.5"),
lty = c(3, 1, 8), lwd = 1.0)
}
powerUtilityPlots()
# Right-Hand Plot Code
quadraticUtilityPlot <- function()
{
v <- seq(0, 1.5, .01)
u <- v - v^2
ylim <- c(-0.7, 0.4)
plot(v, u, type = "l", ylim = ylim, xlab = "v", ylab = "U(v)", lwd = 1.5)
abline(v = .50, lty = "dotted")
abline(h = .25, lty = "dotted")
}
quadraticUtilityPlot()
## Figure 2.35
# Left-Hand Plot Code
rm1 <- 0.18
Beta1 <- c(1.53, 1.36, 1.24, 1.17, 1.06, 0.92, 0.84, 0.76, 0.63, 0.48)
Mu1 <- c(0.26, 0.22, 0.21, 0.21, 0.18, 0.17, 0.16, 0.15, 0.13, 0.12)
Sigma1 <- c(0.49, 0.43, 0.39, 0.37, 0.33, 0.29, 0.27, 0.24, 0.20, 0.17)
SML1 <- rm1 * Beta1
df1 <- data.frame( Beta1, Mu1, Sigma1, SML1)
df1a <- data.frame(x=1, y = rm1)
p1 <- ggplot(df1) +
geom_point(aes(x = Beta1, y = Mu1), color = "black") +
geom_line(aes( x = Beta1, y = SML1), color = "gray20") +
labs(x = "Beta", y = "Mean Excess Return", title = "1931 \u2013 1965") +
theme(plot.title = element_text(hjust = 0.5)) +
annotate(geom="text", x=1.19, y=0.15,
label="CAPM Security Market Line", color="gray20") +
annotate(geom="text", x=.78, y=0.185,
label="Market Portfolio", color="dodgerblue4")
p1 + geom_point(data = df1a, aes(x = x, y = y),
color = "dodgerblue4", shape = 17, size = 3)
# Rigt-Hand Plot Code
rm2 <- 0.08
Beta2 <- c(1.50, 1.30, 1.17, 1.09, 1.03, 0.95, 0.87, 0.78, 0.67, 0.51)
Mu2 <- c(0.06, 0.08, 0.08, 0.08, 0.08, 0.08, 0.07, 0.07, 0.07, 0.06)
Sigma2 <- c(0.31, 0.26, 0.24, 0.22, 0.21, 0.19, 0.18, 0.16, 0.14, 0.12)
SML2 <- rm2 * Beta2
df2 <- data.frame(Beta2, Mu2, Sigma2, SML2)
df2a <- data.frame(x = 1, y = rm2)
p2 <- ggplot(df2) +
geom_point(aes(x = Beta2, y = Mu2), color = "black") +
geom_line(aes( x = Beta2, y = SML2), color = "gray20") +
ylim(0, 0.12) +
labs(x = "Beta", y = "Mean Excess Return", title = "1965 \u2013 1991") +
theme(plot.title = element_text(hjust = 0.5)) +
annotate(geom = "text", x = 1.12, y = 0.055,
label = "CAPM Security Market Line", color = "gray20") +
annotate(geom = "text", x = .80, y = 0.095,
label = "Market Portfolio", color = "dodgerblue4")
p2 + geom_point(data = df2a, aes(x = x, y = y),
color = "dodgerblue4", shape = 17, size = 3)
## Figure 2.36
# Left-Hand Plot Code
ggplot(df1, aes(x = Beta1, y = Sigma1)) +
geom_point() +
geom_smooth(formula = 'y ~ x', method='lm', se = FALSE, linewidth = 0.6, color = "gray20") +
labs(x = "Beta", y = "Standard Deviation") +
geom_text(x = 1.05, y = 0.27, label = "sigma %~~% 0.32 ~ beta", parse=TRUE)
# Right-Hand Plot Code
ggplot(df2, aes(x = Beta2, y = Sigma2)) +
geom_point() +
geom_smooth(formula = 'y ~ x', method='lm', se = FALSE, linewidth = 0.6, color = "gray20") +
labs(x = "Beta", y = "Standard Deviation") +
geom_text(x = 1.05, y = 0.18, label = "sigma %~~% 0.2 ~ beta", parse=TRUE)
## Table 2.7
df <- data.frame(matrix(" ", nrow = 6, ncol = 4))
df$X1 <- c("1/31--12/39", "1/40--12/49", "1/50--12/59", "1/60--12/69", "1/70--12/79","1/80--12/91")
df$X2 <- c(-0.05, 0.03, 0.08, 0.03, 0.01, 0.09)
df$X3 <- c(0.17, 0.10, 0.06, 0.07, 0.10, 0.08)
df$X4 <- c(-0.94, 1.06, 4.25, 1.32, 0.18, 3.90)
#Rename rows and columns and reformat the table
# colnames(df) <- c("Period", "$\\mu_{e}$", "$\\sigma_{e}$", "$t( \\mu )$") # In the text
colnames(df) <- c("Period", "Mean(Excess Return)", "Std. Dev(Excess Return)", "t(Mean(Excess Return))")
df
## Figure 2.37
stocksCRSPweekly <- getPCRAData("stocksCRSPweekly")
dateRange <- c("2004-01-01", "2005-12-31")
stockItems <- c("Date", "TickerLast", "CapGroupLast",
"Return", "MktIndexCRSP", "Ret13WkBill")
returnsAll <- selectCRSPandSPGMI("weekly",
dateRange = dateRange,
stockItems = stockItems,
factorItems = NULL,
subsetType = "CapGroupLast",
subsetValues = "SmallCap",
outputType = "xts")
returns <- returnsAll[ , 1:10]
tsPlotMP(returns, scaleType = "free",layout = c(2,5),
stripText.cex = .45, axis.cex = 0.4,lwd = 0.5)
## Figure 2.38
pspec <- portfolio.spec(assets = names(returns))
pspecFI <- add.constraint(pspec, type = "full_investment")
pspecLO <- add.constraint(portfolio = pspecFI, type = "long_only")
pspecESLO5pct <- add.objective(pspecLO, type = "risk", name = "ES",
arguments = list(p = 0.050))
# Increase n.portfolios below for more accurate vertical dot-dash line
lty <- c("dashed", "solid", "dotted", "dotdash")
col <- c("red", "black", "darkgreen", "darkgreen")
chart.EfficientFrontierCompare(returns, pspecESLO5pct, risk_type = "ES",
guideline = TRUE, cex.axis = 1.2,
match.col = c("StdDev", "ES"),
n.portfolios = 10,
lwd = c(1.3, 1.4, 1.3, 1.0),
col = col, lty = lty,
xlim = c(0.02, 0.08), ylim = c(0.0, 0.012),
legend.loc = "topleft", main = NULL)
## Figure 2.39
chart.EfficientFrontierCompare(returns, pspecESLO5pct, risk_type = "StdDev",
guideline = TRUE, cex.axis = 1.2,
match.col = c("ES", "StdDev"),
n.portfolios = 10,
lwd = c(1.3, 1.4, 1.3, 1.0),
col = col, lty = lty,
xlim = c(0.01, 0.06), ylim = c(0.0, 0.012),
legend.loc = "topleft", main = NULL)
## Figure 2.40
x <- seq(-4.9, 4.9, by = 0.001)
ccopt <- computeTuningPsi_modOpt(0.95)
plot(x, wgt_modOpt(x, cc = ccopt), type = "l",
xlab = "x", ylab = "", cex = 1.5, cex.lab = 1.5)
## Figure 2.41
data(edhec)
hfnames <- c("CA", "CTA", "DIS", "EM", "EMN", "ED", "FIA",
"GM", "LSE", "MA", "RV", "SS", "FOF")
names(edhec) <- hfnames
retLongFIA <- edhec[, "FIA"]
retFIA <- retLongFIA['1998-01-31/1999-12-31', ]
index(retFIA) <- as.yearmon(index(retFIA))
# Plot FIA returns with sample mean and robust mean
mu <- 100*mean(retFIA)
se.mu <- 100*sd(retFIA)/sqrt(24)
x <- RobStatTM::locScaleM(retFIA, eff = .95)
muRob <- 100*x$mu
se.muRob <- 100*x$std.mu
plot.zoo(retFIA, type ="b", xlab = "", ylab = "FIA Returns")
abline(h = muRob/100, col = "blue")
abline(h = mu/100, lty = "dashed", col ="red")
legend(1999.2, -.03, legend = c("Robust Mean", "Sample Mean"), lwd = c(1, 2),
lty = c("solid", "dashed"), col = c("blue", "red"), bty = "n", cex = 1.3)
## Table 2.8
tstat.mu <- mu/se.mu
tstat.muRob <- muRob/se.muRob
SR.classic <- tstat.mu/sqrt(24)
SR.Rob <- tstat.muRob/sqrt(24)
row1 <- round(c(mu, se.mu, tstat.mu, SR.classic), 2)
row2 <- round(c(muRob, se.muRob, tstat.muRob, SR.Rob), 2)
meanEsts <- data.frame(rbind(row1, row2))
names(meanEsts) <- c("Estimate (%)", "Std. Error (%)", "t-Stat", "Sharpe Ratio")
row.names(meanEsts) <- c("Sample Mean", "Robust Mean")
meanEsts
## Figure 2.42
# mOpt 95% Efficiency Tuning Constant
ccModOpt <- computeTuningPsi_modOpt(0.95)
# mOpt M-scale weight function
wgt_mOptScale <- function(x) {
rho_modOpt(x, cc = ccModOpt)/(x^2)}
wgt_mOptScaleMax <- wgt_mOptScale(0.0001)
# Plot mOptScale Weight Function
x <- seq(-5.5, 5.5, 0.01)
ylim <- c(0, 1.4)
ylab <- "w_scale,mOpt(x)"
plot(x,wgt_mOptScale(x)/wgt_mOptScaleMax, ylim = ylim,
ylab = ylab, type = "l", cex.lab = 1.1)
abline(h = 1.0, lty = "dotted")
## Figure 2.43
data(edhec, package = "PerformanceAnalytics")
hfnames <- c("CA", "CTA", "DIS", "EM", "EMN", "ED", "FIA",
"GM", "LSE", "MA", "RV", "SS", "FOF")
names(edhec) <- hfnames
range(index(edhec))
edhec <- edhec[ , 1:12]
returns <- 100*edhec['1998-01-31/1999-12-31']
# Plot edhec returns for 1998-1999
tsPlotMP(returns, type = "l", stripText.cex = 0.7, axis.cex = 0.7)
## Table 2.9
StdDev <- apply(returns,2,sd)
MADM <- apply(returns,2,mad)
resid <- returns - median(returns)
RobSD <- apply(resid, 2, scaleM, family = "mopt")
SDestsMat <- cbind(StdDev, MADM, RobSD)
SDestsMat <- round(SDestsMat,2)
SDests <- data.frame(hfnames[1:12], SDestsMat)
names(SDests) <- c("HFindex", "StdDev", "MADM", "RobSD")
row.names(SDests) <- NULL
SDests
Try the PCRA package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.