Skip to content
Commits on Source (3)
Show whitespace changes
Inline Side-by-side
DESCRIPTION
View file @ 76c90b28
Package: CLA
Version: 0.90-1
Date: 2018-01-25
Date: 2018-02-05
Title: Critical Line Algorithm in Pure R
Author: Yanhao Shi <syhelena@163.com>,
Martin Maechler <maechler@stat.math.ethz.ch>
......
R/findSigMu.R
View file @ 76c90b28
## From Yanhao Shi's thesis -- R/Functions/CLA_analysis.R
#### From Yanhao Shi's thesis -- R/Functions/CLA_analysis.R
findSig <- function(Mu0, result, covar, equal.tol){ # equal.tol > 1e-16
## Vectorized in first argument by Martin Maechler
### Version 1 -- of the code, simply "vectorize via lapply" -- not so fast
findMu <- function(Sig0, result, covar, tol.unir = 1e-6, equal.tol = 1e-6) {
R <- lapply(Sig0, findMu.1,
result=result, covar=covar, tol.unir=tol.unir, equal.tol=equal.tol)
## each R[[i]] has (Mu, weight)
list(Mu = vapply(R, `[[`, numeric(1), "Mu"),
weight = vapply(R, `[[`, numeric(nrow(covar)), "weight"))
}
findSig <- function(Mu0, result, covar, equal.tol = 1e-6) { # equal.tol > 1e-16
R <- lapply(Mu0, findSig.1,
result=result, covar=covar, equal.tol=equal.tol)
## each R[[i]] has (Sig, weight):
list(Sig = vapply(R, `[[`, numeric(1), "Sig"),
weight = vapply(R, `[[`, numeric(nrow(covar)), "weight"))
}
##' Now simplified and very close to findSig0() -- just returning weights additionally
findSig.1 <- function(Mu0, result, covar, equal.tol) { # equal.tol > 1e-16
stopifnot(length(Mu0) == 1L)
ms.w <- result$MS_weight
nd <- order(ms.w[, "Mu"])
n <- length(nd)
## revert order for all (w, mu, sig):
ms.w <- ms.w[nd, ]
w <- result$weights_set[, nd]
ind.uni <- c(TRUE, sapply(1:(n-1), function(i) !
isTRUE(all.equal(ms.w[i, ], ms.w[i+1, ], tol = equal.tol))))
p <- rep(1, n) # pointer, different weights of same ms have same pointer
k <- 1
for(i in 1:n){
if(ind.uni[i]){
p[i] <- k
k <- k + 1
}else{
p[i] <- p[i-1]
}
}
mu.w <- ms.w[ind.uni, "Mu"]
sig.w <- ms.w[ind.uni, "Sig"]
if(Mu0 < min(mu.w) || Mu0 > max(mu.w)) stop(sprintf("Mu0 must be in [%g, %g]",
max(mu.w), min(mu.w)))
mu.w <- ms.w[, "Mu"]
sig.w <- ms.w[, "Sig"]
if(Mu0 < min(mu.w) || Mu0 > max(mu.w))
stop(sprintf("Mu0 must be in [%g, %g]", min(mu.w), max(mu.w)))
i <- findInterval(Mu0, mu.w) # [..)...[..]
if(i == n) { # Mu0 is max(mu.w)
if(i == n || # Mu0 is max(mu.w)
isTRUE(all.equal(mu.w[i], mu.w[i+1], tol = equal.tol))) {
list(Sig = sig.w[i], weight = w[, i])
} else {
m1 <- which(p == i)
m2 <- which(p == i+1)
a <- (Mu0 - mu.w[i+1])/(mu.w[i] - mu.w[i+1])
# solve for a : Mu0 = a* Mu1 + (1-a)* Mu2 :
wm <- matrix(0, n, m1*m2)
for(s in m1){
wm[,(s-1)*m2 + 1:m2] <- a*w[, s] + (1-a)*w[, m2]
}
list(Sig = sqrt(colSums(wm *(covar %*% wm))), #return Sig0
weight = wm) # might have replicate columns
wm <- a*w[, i] + (1-a)*w[, i+1]
list(Sig = sqrt(colSums(wm *(covar %*% wm) )),
weight = wm)
}
}
##' This is much cleaner and logical than the original findSig() ==> use a version of it
##' A version of findSig() only used inside findMu():
findSig0 <- function(Mu0, result, covar, equal.tol = 1e-8){ # equal.tol > 1e-16
stopifnot(length(Mu0) == 1L)
ms.w <- result$MS_weight
nd <- order(ms.w[, "Mu"])
n <- length(nd)
## revert order for all (w, mu, sig):
ms.w <- ms.w[nd, ]
w <- result$weights_set[, nd]
mu.w <- ms.w[, "Mu"]
sig.w <- ms.w[, "Sig"]
i <- findInterval(Mu0, mu.w) # [..)...[..]
if(i == n){ # Mu0 is max(mu.w)
sig.w[i]
} else if(isTRUE(all.equal(mu.w[i], mu.w[i+1], tol = equal.tol))){
if(i == n || isTRUE(all.equal(mu.w[i], mu.w[i+1], tol = equal.tol))) {
sig.w[i]
} else {
a <- (Mu0 - mu.w[i+1])/(mu.w[i] - mu.w[i+1])
......@@ -69,21 +70,21 @@ findSig0 <- function(Mu0, result, covar, equal.tol = 1e-8){ # equal.tol > 1e-16
}
}
findMu <- function(Sig0, result, covar, tol.unir = 1e-6, equal.tol = 1e-6){
findMu.1 <- function(Sig0, result, covar, tol.unir = 1e-6, equal.tol = 1e-6) {
stopifnot(length(Sig0) == 1L)
ms.w <- result$MS_weight
nd <- order(ms.w[, "Sig"])
n <- length(nd)
mu.w <- ms.w[nd, "Mu"]
sig.w <- ms.w[nd, "Sig"]
w <- result$weights_set[, nd]
if(Sig0 < min(sig.w) || Sig0 > max(sig.w)) stop(sprintf("Sig0 must be in [%g, %g]",
max(sig.w), min(sig.w)))
if(Sig0 < min(sig.w) || Sig0 > max(sig.w))
stop(sprintf("Sig0 must be in [%g, %g]", min(sig.w), max(sig.w)))
i <- findInterval(Sig0, sig.w)
if(i == n){ # Mu0 is max(mu.w)
if(i == n || # Mu0 is max(mu.w)
isTRUE(all.equal(mu.w[i], mu.w[i+1], tol = equal.tol))) {
list(Mu = mu.w[i], weight = w[, i])
}else if(isTRUE(all.equal(mu.w[i], mu.w[i+1], tol = equal.tol))){
list(Mu = mu.w[i], weight = w[, i])
} else{
} else { ## FIXME: here are using default equal.tol = 1e-8 in findSig0() !
r <- uniroot(function(mu) findSig0(mu, result, covar) - Sig0,
interval = mu.w[c(i, i+1)], tol=tol.unir)
Mu0 <- r$root
......@@ -94,6 +95,54 @@ findMu <- function(Sig0, result, covar, tol.unir = 1e-6, equal.tol = 1e-6){
}
}
### Version 2 --- vectorize inside -- only those parts that need it
## not yet !!! === For first part, look at
## =========== findSigMu.R+
## ============
if(FALSE) ## will be .... not yet
findMu <- function(Sig0, result, covar, tol.unir = 1e-6, equal.tol = 1e-6){
ms.w <- result$MS_weight
nd <- order(ms.w[, "Sig"])
n <- length(nd)
mu.w <- ms.w[nd, "Mu"]
sig.w <- ms.w[nd, "Sig"]
w <- result$weights_set[, nd]
if(Sig0 < min(sig.w) || Sig0 > max(sig.w))
stop(sprintf("Sig0 must be in [%g, %g]", min(sig.w), max(sig.w)))
ini <- findInterval(Sig0, sig.w)
m <- length(Sig0)
mu0 <- numeric(m)
wt0 <- matrix(NA_real_, n, m)
iBnd <- vapply(ini, function(i) {
(i == n || # Mu0 is max(mu.w)
isTRUE(all.equal(mu.w[i], mu.w[i+1], tol = equal.tol))) # duplicate turning pt
}, NA)
if(any(iBnd)) {
i <- ini[iBnd]
mu0[ iBnd] <- mu.w[i]
wt0[,iBnd] <- w[, i]
}
if(any(iIn <- !iBnd)) { # regular case
i <- ini[iIn]
mus <- vapply(..., function(i) {
r <- uniroot(function(mu) findSig0(mu, result, covar) - Sig0[i],
interval = mu.w[c(i, i+1)], tol = tol.unir)
## TODO check convergence?
r$root
}, numeric(1))
## solve for a : mu = a* Mu1 + (1-a)* Mu2 ==>
a <- (Mu0 - mu.w[i+1])/(mu.w[i] - mu.w[i+1])
w0 <- a* w[, i] + (1-a)*w[, i+1]
mu0[ii] <- mu.w[i]
wt0[,ii] <- w[, i]
}
list(Mu = Mu0, weight = w0)
}
......
TODO
View file @ 76c90b28
......@@ -2,11 +2,11 @@
* ASAP (no longer "Before release of package")
** check arguments e.g., lB <= uB, sum upper Bounds >= 1
** TODO References --> (partly done)
** TODO References --> (mostly done)
*** DONE 1) References from the thesis, including the "buried" python-paper with *WRONG* algo
*** TODO 2) Master thesis: I'd like the thesis to be on our web page
** DONE CLA() should return a (S3) class, "CLA" w/ print() and plot() methods
** TODO findMu() and findSig() regression check examples *before* changing R/findSigMu.R
** TODO findMu() and findSig() regression check examples *before* much changing R/findSigMu.R
** TODO Improve plot() method, using hyperbolic interpolation see R/CLA.R man/plot.CLA.Rd
** TODO A. Norring's Masters thesis has a small 12-asset example (from a published source).
We should add that as a minimally small data set to use in examples,
......
man/CLA.Rd
View file @ 76c90b28
......@@ -110,4 +110,6 @@ if(require(Matrix)) { ## visualize how weights change "along turning points"
xlab = "turning point", ylab = "asset number"))
}
}
\keyword{optimize}
\keyword{arith}
man/findMu.Rd
View file @ 76c90b28
\name{findMu}
\alias{findMu} %% keep in sync with ./findSig.Rd
%% ^^^^^^^^^^ %% ??? TODO ??? one help page also for findSig() ???
%% ^^^^^^^^^^ ?? TODO ?? one help page also for findSig() ???
\title{Find mu(W) and W, given sigma(W) and CLA result}
\description{
Find \eqn{\mu(W)} and \eqn{W}, given \eqn{\sigma(W}) and \code{\link{CLA}} result.
Find \eqn{\mu(W)} and \eqn{W}, given \eqn{\sigma(W}) and
\code{\link{CLA}} result.
}
\usage{
findMu(Sig0, result, covar, tol.unir = 1e-06, equal.tol = 1e-06)
}
\arguments{
\item{Sig0}{numeric vector of \eqn{\sigma(W)} values.}
\item{result}{ .. }
\item{covar}{ ... }
\item{result}{a \code{\link{list}} with components \code{MS_weight}
and \code{weights_set} as resulting from \code{\link{CLA}()}.}
\item{covar}{the same \eqn{n \times n}{n x n} covariance matrix (of
asset returns) as the argument of \code{\link{CLA}()}.}
\item{tol.unir}{numeric tolerance passed to \code{\link{uniroot}}.}
\item{equal.tol}{numeric tolerance to be used in
\code{\link{all.equal}(.., tolerance = equal.tol)}.}% TODO: is this correct?
\code{\link{all.equal}(.., tolerance = equal.tol)} in the check to
see if the \eqn{\mu}{mu} of two neighbouring turning points are
equal.}
}
\value{
a \code{\link{list}} with components
\item{Mu}{.....}
\item{weight}{.....}
\item{Mu}{numeric vector of same length, say \eqn{M}, as \code{Sig0}.}
\item{weight}{numeric \eqn{n \times M}{n x M} matrix of weights.}
}
\references{
Master thesis, p.33
......@@ -30,7 +35,21 @@ findMu(Sig0, result, covar, tol.unir = 1e-06, equal.tol = 1e-06)
\code{\link{CLA}}, \code{\link{MS}}.
}
\examples{
## TODO
data(muS.sp500)
## Full data taking too much time for example
set.seed(2016)
iS <- sample.int(length(muS.sp500$mu), 17)
cov17 <- muS.sp500$covar[iS, iS]
CLsp.17 <- CLA(muS.sp500$mu[iS], covar=cov17, lB=0, uB = 1/2)
CLsp.17 # 16 turning points
summary(tpS <- CLsp.17$MS_weights[,"Sig"])
str(s0 <- seq(0.0186, 0.0477, by = 0.0001))
mu.. <- findMu(s0, result=CLsp.17, covar=cov17)
str(mu..)
stopifnot(dim(mu..$weight) == c(17, length(s0)))
plot(s0, mu..$Mu, xlab=quote(sigma), ylab = quote(mu),
type = "o", cex = 1/4)
points(CLsp.17$MS_weights, col = "tomato", cex = 1.5)
}
%% \keyword{ ~kwd1 }% use one of RShowDoc("KEYWORDS")
%% \keyword{ ~kwd2 }% __ONLY ONE__ keyword per line
\keyword{optimize}
\keyword{arith}
man/findSig.Rd
View file @ 76c90b28
......@@ -2,37 +2,56 @@
\alias{findSig} %% keep in sync with ./findMu.Rd
\title{Find sigma(W) and W, given mu(W) and CLA result}
\description{
Find \eqn{\sigma(W)} and \eqn{W}, given \eqn{\mu(W}) and \code{\link{CLA}} result.
Find \eqn{\sigma(W)} and \eqn{W}, given \eqn{\mu(W}) and
\code{\link{CLA}} result.
}
\usage{
findSig(Mu0, result, covar, equal.tol)
}
\arguments{
\item{Mu0}{numeric vector}
\item{result}{ .. }
\item{covar}{ ... }
\item{Mu0}{numeric vector of \eqn{\mu(W)} values.}
\item{result}{a \code{\link{list}} with components \code{MS_weight}
and \code{weights_set} as resulting from \code{\link{CLA}()}.}
\item{covar}{the same \eqn{n \times n}{n x n} covariance matrix (of
asset returns) as the argument of \code{\link{CLA}()}.}
\item{equal.tol}{numeric tolerance to be used in
\code{\link{all.equal}(.., tolerance = equal.tol)}.}% TODO: is this correct?
\code{\link{all.equal}(.., tolerance = equal.tol)} in the check to
see if the \eqn{\mu}{mu} of two neighbouring turning points are
equal.}
}
%% \details{
%% %% ~~ If necessary, more details than the description above ~~
%% }
\value{
a \code{\link{list}} with components
\item{Sig}{.....}
\item{weight}{.....}
\item{Sig}{numeric vector of same length, say \eqn{M}, as \code{Mu0}.}
\item{weight}{numeric \eqn{n \times M}{n x M} matrix of weights.}
}
\references{
Master thesis, p.33
}
\seealso{
\code{\link{findMu}},
\code{\link{CLA}}, \code{\link{MS}}.
}
\examples{
## ...
data(muS.sp500)
## Full data taking too much time for example: Subset of n=21:
set.seed(2018)
iS <- sample.int(length(muS.sp500$mu), 21)
cov21 <- muS.sp500$covar[iS, iS]
CLsp.21 <- CLA(muS.sp500$mu[iS], covar=cov21, lB=0, uB = 1/2)
CLsp.21 # 14 turning points
summary(tpM <- CLsp.21$MS_weights[,"Mu"])
str(m0 <- c(min(tpM),seq(0.00205, 0.00525, by = 0.00005), max(tpM)))
sig. <- findSig(m0, result=CLsp.21, covar=cov21)
str(sig.)
stopifnot(dim(sig.$weight) == c(21, length(m0)))
plot(sig.$Sig, m0, xlab=quote(sigma), ylab = quote(mu),
type = "o", cex = 1/4)
points(CLsp.21$MS_weights, col = "tomato", cex = 1.5)
title("Efficient Frontier from CLA()")
mtext("findSig() to interpolate between turning points", side=3)
}
%% \keyword{ ~kwd1 }% use one of RShowDoc("KEYWORDS")
\keyword{arith}
\keyword{dplot}