version 0.90-0

DESCRIPTION

@@ -0,0 +1,20 @@
+Package: CLA
+Version: 0.90-0
+Date: 2017-09-25
+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
+Suggests: fGarch, FRAPO, Matrix
+Description: Implements 'Markovitz' Critical Line Algorithm ('CLA') for
+  classical mean-variance portfolio optimization.  Care has been taken for
+  correctness in light of previous buggy implementations.
+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
+Repository: CRAN
+Date/Publication: 2018-01-15 18:50:03 UTC

MD5

@@ -0,0 +1,22 @@
+13c290742a8eb73c94984c64ab7ca581 *DESCRIPTION
+0f8952946291e6f083553ea5594bf4c4 *LICENSE
+e7ca4aca5e757e488a10dc80aaaacba7 *NAMESPACE
+6807f29369c23e9a7d5d42fb6ae93a17 *R/Auxiliaries.R
+b3f72360e5e381ee50c8cc2de33abe6e *R/CLA.R
+d5f4802f25daad1cc5239a6f69ef6f2d *R/Fenv.R
+5df22ffbfd27cb038c2be63955c4e73d *R/Ver8.R
+40a0acc69c15ad4aeb5cdd1fec3c48d1 *R/Ver9.R
+4a429eee36db470a747697d0cffee992 *R/findSigMu.R
+4dbc9f200fbe9beb38d2c95294aa5395 *R/muSigmaGarch.R
+0e726966c0774a979c68dba1bcc77751 *README.md
+298e672d31caab38d077c0f54f5186e4 *TODO
+c8190d1ec7f8bf4ed04e9dc0357ed50b *data/datalist
+cc397e49de11b7601bfa194e83163d44 *data/muS.sp500.rda
+1a3d508743e85bd627f871a7036b2e6f *man/CLA.Rd
+5c3484d17c7e711b3ae3e7718a22171b *man/MS.Rd
+d96d6da6e1ca33dae117ba62b491cfd8 *man/findMu.Rd
+5ffc02797c2768a2837ce6b1d4723432 *man/findSig.Rd
+58b004ce949485b4713966332cca331f *man/muS.sp500.Rd
+2f1bed365394d2f82ac35e9bf980c3ce *man/muSigmaGarch.Rd
+f04b1c9ad9c9f84b5a753977d8b953be *tests/SP500-ex.R
+dadeb6bbecaab7764ae4bc8a7ce640da *tests/wtsn0.rds

NAMESPACE

@@ -0,0 +1,15 @@
+
+## NOT importing fGarch here -- only depending on using muSigmaGarch()
+
+importFrom("stats",
+           cor, predict, uniroot)
+## importFrom(, ..)
+## not yet on CRAN : importFrom("sfsmisc", funEnv)
+
+export(CLA
+       , MS
+       , findSig, findMu
+       , muSigmaGarch
+       )
+
+S3method(print, CLA)

R/Auxiliaries.R

@@ -0,0 +1,14 @@
+## In next version of  package 'sfsmisc' - for now here (not exported)
+funEnv <- function(..., envir = NULL, parent = parent.frame(),
+                   hash = (...length() > 100), size = max(29L, ...length())) {
+    e <- list2env(list(...), envir=envir, parent=parent, hash=hash, size=size)
+    for(n in names(e)) ## iff function or formula, set environment to 'e':
+	if(is.function(e[[n]]) || (is.call(e[[n]]) &&
+				   inherits(e[[n]], "formula")))
+	    environment(e[[n]]) <- e
+    e
+}
+
+if(!is.function(.BaseNamespaceEnv$...length)) # ...length() only in R >= 3.5.0
+    ## kludgy substitute, using parent.env() -- but it works (sometimes) in funEnv()
+    ...length <- function() eval(quote(length(list(...))), envir = parent.frame())

R/CLA.R

@@ -0,0 +1,197 @@
+#### Started from "Version 8" (Ver8.R)
+
+## Initialize the weight -- Find first free weight
+initAlgo <- function(mu, lB, uB ) {
+    ## New-ordered return, lB, uB with decreasing return
+    w <- c()
+    index.new <- order(mu,decreasing = TRUE) # new order with decreasing return
+    lB.new <- lB[index.new]
+    uB.new <- uB[index.new]
+    ## free weight - starting solution
+    i.new <- 0
+    w.new <- lB.new # initially
+    while(sum(w.new) < 1) {
+        i.new <- i.new + 1
+        w.new[i.new] <- uB.new[i.new]
+    }
+    w.new[i.new] <- 1 - sum(w.new[-i.new])
+    w[index.new] <- w.new                # back to original order
+    i <- index.new[i.new]
+    ## return the index of first free asset and vector w :
+    list(index = i, weights = w)
+}
+
+## getMatrices -----------------------------------------------------------------
+getMatrices <- function(mu, covar, w, f) {
+    ## Slice covarF,covarFB,covarB,muF,muB,wF,wB
+    covarF <- covar[f,f]
+    muF <- mu[f]
+    b <- seq_along(mu)[-f]
+    covarFB <- covar[f,b]
+    wB <- w[b]
+    list(covarF = covarF, covarFB = covarFB, muF = muF, wB = wB)
+}
+
+computeInv <- function(get) {
+    solve(get$covarF, cbind(1, get$muF, get$covarFB %*% get$wB, deparse.level = 0L))
+}
+
+## computeW -----------------------------------------------------------------
+computeW <- function(lam, inv, wB) {
+    ## w2 <- inv[,1]; w3 <- inv[,2]; w1 <- inv[,3]
+    inv.s <- colSums(inv) # g1 <- inv.s[2]; g2 <- inv.s[1]; g4 <- inv.s[3]
+    ## 1) compute gamma
+    g <- (-lam * inv.s[2] + (1- sum(wB) + inv.s[3]))/inv.s[1]
+    ## 2) compute free weights
+    list(wF = - inv[,3] + g * inv[,1] + lam * inv[,2], gamma = g)
+}
+
+## computeLambda --------------------------------------------------------------
+computeLambda <- function(wB, inv, i, bi.input) {
+    inv.s <- colSums(inv)
+    ## c1 <- inv.s[1]; l2 <- inv.s[3]; c2i <- inv[i,2];
+    ## c3 <- inv.s[2]; c4i <- inv[i, 1]; l1 <- sum(wB)
+    c1 <- inv.s[1]
+    if(length(bi.input) == 1L) {                    # 1.bound to free
+        c4i <- inv[i, 1]
+        Ci <- - c1 * inv[i, 2] + inv.s[2] * c4i
+        ## return lambda :
+        if(Ci == 0)
+            0
+        else
+            ((1- sum(wB) + inv.s[3])* c4i- c1 * (bi.input + inv[i, 3]))/Ci
+    }
+    else {                                          # 2.free to bound
+        c4i <- inv[, 1]
+        Ci <- - c1 * inv[, 2] + inv.s[2] * c4i
+        bi          <- bi.input[i, 1]         # bi.lB
+        bi[Ci > 0]  <- bi.input[i[Ci > 0], 2] # bi.uB
+        bi[Ci == 0] <- 0
+        ## return lambda and boundary :
+        list(lambda = ((1- sum(wB) + inv.s[3]) * c4i- c1 *(bi + inv[, 3]))/Ci,
+             bi = bi)
+    }
+}
+
+MS <- function(weights_set, mu, covar) {
+    Sig2 <- colSums(weights_set *(covar %*% weights_set) )
+    cbind(Sig = sqrt(Sig2), Mu = as.vector(t(weights_set) %*% mu))
+}
+
+
+CLA <- function(mu, covar, lB, uB, tol.lambda = 1e-7,
+                give.MS = TRUE, keep.names = TRUE) {
+    ## minimal argument checks
+    cl <- match.call()
+    n <- length(mu)
+    if(length(lB) == 1) lB <- rep.int(lB, n)
+    if(length(uB) == 1) uB <- rep.int(uB, n)
+    stopifnot(is.numeric(mu), is.matrix(covar), identical(dim(covar), c(n,n)),
+              is.numeric(lB), length(lB) == n,
+              is.numeric(uB), length(uB) == n, lB <= uB)# and in [0,1]
+    ## Compute the turning points, free sets and weights
+    ans <- initAlgo(mu, lB, uB)
+    f <- ans$index
+    w <- ans$weights
+    weights_set <- w  # store solution
+    lambdas <- NA  # The first step has no lambda or gamma, add NA instead.
+    gammas <- NA
+    free_indices <- list(f)
+    lam <- 1 # set non-zero lam
+    while ( lam > 0 && length(f) < length(mu)) {
+      ## 1) case a): Bound one free weight  F -> B
+      l_in <- 0
+      if(length(f) > 1 ) {
+        compl <- computeLambda(wB = w[-f], inv = inv, # inv from last step k (k > 1)
+                               i = f, bi.input = cbind(lB, uB))
+        lam_in <- compl$lambda
+        bi     <- compl$bi
+        k <- which.max(lam_in)
+        i_in  <-      f[k]
+        bi_in <-     bi[k]
+        l_in  <- lam_in[k]
+      }
+
+      ## 2) case b): Free one bounded weight  B -> F
+      b <- seq_along(mu)[-f]
+      inv_list <- lapply(b, function(bi) {
+        get_i <- getMatrices(mu, covar, w, c(f,bi))
+        computeInv(get_i)
+      })
+
+      fi <- length(f) + 1
+      lam_out <- sapply(seq_along(b), function(i) {
+        computeLambda(wB = w[b[-i]], inv = inv_list[[i]],
+                      i = fi, bi.input = w[b[i]])
+      })
+
+      if (length(lambdas) > 1 && any(!(sml <- lam_out < lam*(1-tol.lambda)))) {
+        lam_out <- lam_out[sml]
+        b       <- b      [sml]
+        inv_list <- inv_list[sml]
+      }
+      k <- which.max(lam_out)
+      i_out <- b      [k] # one only !
+      l_out <- lam_out[k]
+      inv_out <- inv_list[[k]]
+
+      ## 3) decide lambda
+      lam <- max(l_in, l_out, 0)
+      if(lam > 0) { # remove i_in from f; or add i_out into f
+        if(l_in > l_out ) {
+          w[i_in] <- bi_in  # set value at the correct boundary
+          f <- f[f != i_in]
+          getM <- getMatrices(mu, covar, w, f)
+          inv <- computeInv(getM)
+        }
+        else {
+          f <- c(f,i_out)
+          inv <- inv_out
+        }
+        compW <- computeW(lam, inv = inv, wB = w[-f])
+      }
+      else{ # 4) if max(l_in, l_out) < 0, "stop" when at the min var solution!
+        compW <- computeW(lam = lam, inv = inv, wB = w[-f])
+        ## muF = 0 not necessary, get1 replaced by getM (ie getM from previous step)
+      }
+
+      wF <- compW$wF
+      g <- compW$gamma
+      w[f] <- wF[seq_along(f)]
+
+      lambdas <- c(lambdas, lam)
+      weights_set <- cbind(weights_set, w, deparse.level = 0L) # store solution
+      gammas <- c(gammas, g)
+      free_indices <- c(free_indices, list(sort(f)))
+    }# end While
+
+    if(keep.names) rownames(weights_set) <- names(mu)
+#### FIXME: This change also needs changes in the ../tests/SP500-ex.R !
+    structure(class = "CLA",
+              list(weights_set = weights_set,
+                   free_indices = free_indices,
+                   gammas = gammas, lambdas = lambdas,
+                   MS_weights = if(give.MS)
+                           MS(weights_set = weights_set, mu = mu, covar = covar),
+                   call = cl))
+}
+
+
+## print method
+print.CLA <- function(x, ...) {
+    cat("Call: ", paste(deparse(x$call), sep = "\n", collapse = "\n"),
+        "\n", sep = "")
+    wts <- x$weights_set
+    n <- nrow(wts)
+    nT <- ncol(wts)
+    cat(gettextf("Critical Line Algorithm for n = %d assets, ", n),
+        gettextf("resulting in  %d  turning points", nT),"\n", sep="")
+    ## For now; we can do better later:
+    cat("Overview of result parts:\n")
+    utils::str(x[1:5], max.level = 1, give.attr=FALSE)
+    ## TODO better, e.g., summarizing the number "active assets"
+    ## and those with non-0 weights -- only if lower bounds were (mostly) 0
+    invisible(x)
+}
+
+## TODO:  plot method -- efficient frontier