Commit d5fedbe9 by Martin Maechler

add examples to find*(); also \keyword{}s

parent 759f86e5
 ## 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")) } ## This does _not_ work (at least in some cases) 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) ... ... @@ -25,7 +48,7 @@ findSig <- function(Mu0, result, covar, equal.tol){ # equal.tol > 1e-16 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))) min(mu.w), max(mu.w))) i <- findInterval(Mu0, mu.w) # [..)...[..] if(i == n) { # Mu0 is max(mu.w) list(Sig = sig.w[i], weight = w[, i]) ... ... @@ -35,8 +58,8 @@ findSig <- function(Mu0, result, covar, equal.tol){ # equal.tol > 1e-16 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] 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 ... ... @@ -44,8 +67,10 @@ findSig <- function(Mu0, result, covar, equal.tol){ # equal.tol > 1e-16 } ##' This is much cleaner and logical than findSig.1 above ==> use it everywhere ##' 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) ... ... @@ -57,11 +82,10 @@ findSig0 <- function(Mu0, result, covar, equal.tol = 1e-8){ # equal.tol > 1e-16 mu.w <- ms.w[ , "Mu"] sig.w <- ms.w[ , "Sig"] 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))) { sig.w[i] } else if(isTRUE(all.equal(mu.w[i], mu.w[i+1], tol = equal.tol))){ sig.w[i] }else{ } else { a <- (Mu0 - mu.w[i+1])/(mu.w[i] - mu.w[i+1]) # solve for a : Mu0 = a* Mu1 + (1-a)* Mu2 : w0 <- a*w[, i] + (1-a)*w[, i+1] ... ... @@ -69,21 +93,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"] 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 +118,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) } ... ...
 ... ... @@ -110,4 +110,6 @@ if(require(Matrix)) { ## visualize how weights change "along turning points" xlab = "turning point", ylab = "asset number")) } } \keyword{optimize} \keyword{arith}
 \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}  ... ... @@ -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)) 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}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!