version 0.90-1

cran · Feb 6, 2018 · cb3db2a · cb3db2a
1 parent 06db33b
commit cb3db2a
Show file tree

Hide file tree

Showing 13 changed files with 381 additions and 113 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,12 +1,12 @@
 Package: CLA
-Version: 0.90-0
-Date: 2017-09-25
+Version: 0.90-1
+Date: 2018-02-05
 Title: Critical Line Algorithm in Pure R
 Author: Yanhao Shi <syhelena@163.com>,
 	Martin Maechler <maechler@stat.math.ethz.ch>
 Maintainer: Martin Maechler <maechler@stat.math.ethz.ch>
 Depends: R (>= 3.2.0)
-Imports: stats, graphics
+Imports: stats, grDevices, graphics
 Suggests: fGarch, FRAPO, Matrix
 Description: Implements 'Markovitz' Critical Line Algorithm ('CLA') for
   classical mean-variance portfolio optimization.  Care has been taken for
@@ -15,6 +15,6 @@ License: GPL (>= 3) | file LICENSE
 Encoding: UTF-8
 URL: https://gitlab.math.ethz.ch/maechler/CLA/
 NeedsCompilation: no
-Packaged: 2018-01-12 17:29:26 UTC; maechler
+Packaged: 2018-02-05 17:10:03 UTC; maechler
 Repository: CRAN
-Date/Publication: 2018-01-15 18:50:03 UTC
+Date/Publication: 2018-02-06 09:51:58 UTC
diff --git a/MD5 b/MD5
@@ -1,22 +1,24 @@
-13c290742a8eb73c94984c64ab7ca581 *DESCRIPTION
+512a7cee29c2c62071562ba8475b2058 *DESCRIPTION
 0f8952946291e6f083553ea5594bf4c4 *LICENSE
-e7ca4aca5e757e488a10dc80aaaacba7 *NAMESPACE
+22fd307566a5e9439cce1f4c0291d23f *NAMESPACE
 6807f29369c23e9a7d5d42fb6ae93a17 *R/Auxiliaries.R
-b3f72360e5e381ee50c8cc2de33abe6e *R/CLA.R
+593a64eb2829e6e9f6fb008de7a75569 *R/CLA.R
 d5f4802f25daad1cc5239a6f69ef6f2d *R/Fenv.R
 5df22ffbfd27cb038c2be63955c4e73d *R/Ver8.R
 40a0acc69c15ad4aeb5cdd1fec3c48d1 *R/Ver9.R
-4a429eee36db470a747697d0cffee992 *R/findSigMu.R
+70cac219599038fb4f2fe1fc74526d76 *R/findSigMu.R
 4dbc9f200fbe9beb38d2c95294aa5395 *R/muSigmaGarch.R
 0e726966c0774a979c68dba1bcc77751 *README.md
-298e672d31caab38d077c0f54f5186e4 *TODO
+d4511d32e82f69e3e68db03aa2352394 *TODO
+dad7eb75ba23663ef34faf2a7a2bd362 *build/partial.rdb
 c8190d1ec7f8bf4ed04e9dc0357ed50b *data/datalist
 cc397e49de11b7601bfa194e83163d44 *data/muS.sp500.rda
-1a3d508743e85bd627f871a7036b2e6f *man/CLA.Rd
-5c3484d17c7e711b3ae3e7718a22171b *man/MS.Rd
-d96d6da6e1ca33dae117ba62b491cfd8 *man/findMu.Rd
-5ffc02797c2768a2837ce6b1d4723432 *man/findSig.Rd
+a33eb52df743009652e8d924272e8ac7 *man/CLA.Rd
+280c324facad3d68412e7897355e94d8 *man/MS.Rd
+f559bc67954fbb61b0e1ba92e699a945 *man/findMu.Rd
+e69d79271673c26c3297ba705d9ac69e *man/findSig.Rd
 58b004ce949485b4713966332cca331f *man/muS.sp500.Rd
 2f1bed365394d2f82ac35e9bf980c3ce *man/muSigmaGarch.Rd
-f04b1c9ad9c9f84b5a753977d8b953be *tests/SP500-ex.R
+aa0f96e0539524b42905e1eab6bb0939 *man/plot.CLA.Rd
+036102fe91522915c6e10bffb8b77c4a *tests/SP500-ex.R
 dadeb6bbecaab7764ae4bc8a7ce640da *tests/wtsn0.rds
diff --git a/NAMESPACE b/NAMESPACE
@@ -3,13 +3,19 @@
 
 importFrom("stats",
            cor, predict, uniroot)
-## importFrom(, ..)
+
+importFrom("grDevices", adjustcolor)
+
+importFrom("graphics", plot)# including the plot.default() method
+
 ## not yet on CRAN : importFrom("sfsmisc", funEnv)
 
+## -------------- CLA  Exports -----------------------
 export(CLA
        , MS
        , findSig, findMu
        , muSigmaGarch
        )
 
 S3method(print, CLA)
+S3method(plot, CLA)
diff --git a/R/CLA.R b/R/CLA.R
@@ -194,4 +194,36 @@ print.CLA <- function(x, ...) {
     invisible(x)
 }
 
-## TODO:  plot method -- efficient frontier
+### TODO:  plot method -- efficient frontier
+
+## As basically from  .../YanhaoShi/R/Functions/Plot.R :
+MS_plot <- function(ms, type = "o",
+                    main = "Efficient Frontier",
+                    xlab = expression(sigma(w)),
+                    ylab = expression(mu(w)),
+                    col = adjustcolor("blue", alpha.f = 0.5),
+                    pch = 16, ...) {
+    ## list of weights_set, legend...
+    stopifnot(is.matrix(ms), ncol(ms) == 2)
+    plot(ms[,"Sig"], ms[,"Mu"], type=type, pch=pch, col=col,
+         xlab = xlab, ylab=ylab, main=main, ...)
+}
+
+
+## FIXME: --> see also in ../man/plot.CLA.Rd
+## -----
+## 0) Use findMu() and findSig() to draw the lines *between*
+## 1) Learn from Tobias Setz to plot the lower part of the feasible region
+## 2) Better title, using 'call'
+## 3) mark some critical points particularly
+## 4) give information about the *number* critical points / weights sets
+## 5) consider using a  'add = FALSE' argument and then use 'lines()'
+plot.CLA <- function(x, type = "o", main = "Efficient Frontier",
+                    xlab = expression(sigma(w)),
+                    ylab = expression(mu(w)),
+                    col = adjustcolor("blue", alpha.f = 0.5),
+                    pch = 16, ...) {
+    stopifnot(is.matrix(ms <- x$MS_weights))
+    plot(ms[,"Sig"], ms[,"Mu"], type=type, pch=pch, col=col,
+         xlab=xlab, ylab=ylab, main=main, ...)
+}
diff --git a/R/findSigMu.R b/R/findSigMu.R
@@ -1,89 +1,90 @@
-## 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)
+  } else {
     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"]
+  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{
+  } 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]
     sqrt(colSums(w0 *(covar %*% w0) )) #return Sig0
   }
 }
 
-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)
-    list(Mu = mu.w[i], weight = w[, i])
-  }else if(isTRUE(all.equal(mu.w[i], mu.w[i+1], tol = equal.tol))){
+  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{
+  } 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)
+}
+
 
 
 

diff --git a/TODO b/TODO
@@ -1,6 +1,16 @@
 ##-*- org -*-->  Emacs .. [Tab] key + [Org] menu; C-c C-o to follow links
 
-* Before release of package
-** DONE CLA() should return a (S3) class, "CLA"
-*** TODO --> print() and plot() (S3) methods**** plot(): plot efficient frontier
+* ASAP (no longer "Before release of package")
 ** check arguments  e.g., lB <= uB, sum upper Bounds >= 1
+** 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* 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,
+   e.g. plot(). His thesis is in ~/Betreute-Arbeiten/YanhaoShi/Previous_Work/
+
+* With more time, also, e.g., for a short R Journal paper
+** SparseMatrix plot of the weights
diff --git a/build/partial.rdb b/build/partial.rdb