MPT

La Théorie Moderne du Portefeuille (MPT) de H. Markowitz en quelques lignes de code.


rm(list=ls())

# -------------------------------------------------------------------------- #
# yHist
# Pour télécharger des données sur ichart.yahoo.com
# -------------------------------------------------------------------------- #

yHist = function(sym, from, to = Sys.Date(), freq = "d") {
      f <- as.Date(from[1])
      t <- as.Date(to[1])
      q <- tolower(substr(freq[1], 1, 1))
      if(! q %in% c("d", "w", "m")) stop("! freq %in% c('d', 'w', 'm')")
      cn <- c("Sym", "Date", "Adj.Close")
      if(length(sym) > 1) {
            k <- lapply(sym, yHist, from = from, to = to, freq = freq)
            res <- do.call(rbind, k)
      } else {
            url <- paste(c(
                 "http://ichart.yahoo.com/table.csv?",
                  "s=", sym,
                  "&a=", as.integer(format(f, "%m")) - 1,
                  "&b=", as.integer(format(f, "%d")),
                  "&c=", as.integer(format(f, "%Y")),
                  "&d=", as.integer(format(t, "%m")) - 1,
                  "&e=", as.integer(format(t, "%d")),
                  "&f=", as.integer(format(t, "%Y")),
                  "&g=", q,
                  "&ignore=.csv"), collapse = "")
            ans <- try(read.csv(url), silent = TRUE)
            if(class(ans) == "try-error") {
                  res <- NULL
            } else {
                  ans <- ans[order(ans$Date), ]
                  res <- subset(cbind(Sym = sym, ans), select = cn)
                  rownames(res) <- NULL
            }
      }
      return(res)
}

# -------------------------------------------------------------------------- #
# dvar
# Calcule des % de variation.
# -------------------------------------------------------------------------- #

dvar = function(x) {
 if( is.matrix(x) ) {
  res <- apply(x, 2, dvar)
  dimnames(res) <- list(rownames(x)[-1], colnames(x))
 } else {
  res <- x[-1]/x[-length(x)]-1
  names(res) <- names(x)[-1]
 }
 return(res)
}

# ========================================================================== #
# SCRIPT
# ========================================================================== #

# L'indice DJIA (j'ai juste remplacé GS par AAPL).

Sym <- c("AAPL", "AXP", "BA", "CAT", "CSCO", "CVX", "DD", "DIS", "GE", "HD",
 "IBM", "INTC", "JNJ", "JPM", "KO","MCD", "MMM", "MRK", "MSFT", "NKE",
 "PFE", "PG", "T", "TRV", "UNH", "UTX", "V", "VZ", "WMT", "XOM")

# Téléchargement

Dta <- yHist(c(Sym, "^GSPC"), "2013-11-20", "2014-11-20", "d")

# On transforme tout ça en matrice de prix.

P <- with(Dta, tapply(Adj.Close, list(Date, Sym), mean))

# On sauvegarde les dates pour usage ultérieur.

Dates <- as.Date(rownames(P))

# On isole le S&P 500 du reste de nos données.

Ix <- P[, "^GSPC", drop = FALSE]
P <- P[, 1:(ncol(P)-1)]

# Un taux sans risque bricolé rapidement...

Rf <- as.matrix(1.0012^(as.numeric(diff(Dates))/365)-1)

# Matrices de variations quotidiennes

X <- dvar(P)
Xm <- dvar(Ix)

# Portefeuille équipondéré

w <- matrix(1/30, 30, 1)
rownames(w) <- Sym

# Qu'est-ce que ça donne ?
po <- X %*% w

y <- c(1, cumprod(1 + po))
plot(Dates, y, type = "l")

# Voyons la volatilité de ce portefeuille.
sd(po) * sqrt(260)

# ========================================================================== #
# MINIMUM VARIANCE
# ========================================================================== #

# On charge *quadprog*. Si vous ne l'avez pas : install.packages("quadprod")
library("quadprog")

# Matrice de variance-covariance.
S <- cov(X)

# Tout le monde à 0 !
R <- matrix(0, 1, 30)

# Nous voulons que la somme des poids soit égale à 1 ET que tous les poids
# soient supérieurs à 0.
A <- cbind(1, diag(1, 30, 30))
B <- matrix(c(1, rep(0, 30)))

# Attention !
res <- solve.QP(S, R, A, B, 1)

w <- matrix(res$solution, 30, 1)
rownames(w) <- Sym

# 21% de McDo
round(w*100, 2)

po <- X %*% w
sd(po) * sqrt(260)

# Ce portefeuille est imbattable
y <- c(1, cumprod(1 + po))
plot(Dates, y, type = "l")

# ========================================================================== #
# INDICIEL
# ========================================================================== #

Xs <- X - drop(Xm)

# Matrice de variance-covariance.
S <- cov(Xs)

# Tout le monde à 0 !
R <- matrix(0, 1, 30)

# Nous voulons que la somme des poids soit égale à 1 ET que tous les poids
# soient supérieurs à 0.
A <- cbind(1, diag(1, 30, 30))
B <- matrix(c(1, rep(0, 30)))

# Attention !
res <- solve.QP(S, R, A, B, 1)

w <- matrix(res$solution, 30, 1)
rownames(w) <- Sym

po <- X %*% w
sd(po - Xm) * sqrt(260)

# Visualisons
y <- c(1, cumprod(1 + po))
ym <- c(1, cumprod(1 + Xm))
plot(Dates, y, type = "l")
lines(Dates, ym, lty = "dashed")

# ========================================================================== #
# LA FRONTIERE EFFICIENTE
# ========================================================================== #

R <- colMeans(X)
S <- cov(X)
A <- cbind(1, diag(1, 30, 30))
B <- matrix(c(1, rep(0, 30)))

res <- solve.QP(S, R, A, B, 1)

w <- matrix(res$solution, 30, 1)
rownames(w) <- Sym

# LOL ! 100% Apple !
round(w*100, 2)

# Comment faire baisser cette volatilité ?
# Facile : on va utiliser le paramètre q

q <- .5
res <- solve.QP(S, q*R, A, B, 1)

w <- matrix(res$solution, 30, 1)
rownames(w) <- Sym

# Intel rentre dans notre portefeuille.
round(w*100, 2)

po <- X %*% w
sd(po) * sqrt(260)

# Réduisons encore q :

q <- .1
res <- solve.QP(S, q*R, A, B, 1)

w <- matrix(res$solution, 30, 1)
rownames(w) <- Sym

# UnitedHealth rentre dans notre portefeuille.
round(w*100, 2)

po <- X %*% w
sd(po) * sqrt(260)

# Voyons ce que ça donne en faisant varier q entre 0 et .1 :
q <- seq(0, .1, len = 100)

W1 <- matrix(NA, 100, 30)
colnames(W1) <- Sym

Z1 <- matrix(NA, 100, 2)
colnames(Z1) <- c("mean", "sd")

for(i in 1:length(q)) {
 ans <- solve.QP(S, q[i]*R, A, B, 1)
 W1[i, ] <- ans$solution
 po <- X %*% matrix(ans$solution, 30, 1)
 Z1[i, ] <- c(mean(po), sd(po))
}

# La frontière efficiente !
plot(Z1[, "sd"], Z1[, "mean"], type = "b")

# ========================================================================== #
# ACTIF SANS RISQUE
# ========================================================================== #

XX <- cbind(X, Rf)
R <- colMeans(XX)
S <- cov(XX)
A <- cbind(1, diag(1, 31, 31))                 # Attention : 31 !
B <- matrix(c(1, rep(0, 31)))                  # Attention : 31 !

q <- seq(0, .1, len = 100)

W2 <- matrix(NA, 100, 31)                      # Attention : 31 !
colnames(W2) <- Sym

Z2 <- matrix(NA, 100, 2)
colnames(Z2) <- c("mean", "sd")

for(i in 1:length(q)) {
 ans <- solve.QP(S, q[i]*R, A, B, 1)
 W2[i, ] <- ans$solution
 po <- XX %*% matrix(ans$solution, 31, 1) # Attention : 31 !
 Z2[i, ] <- c(mean(po), sd(po))
}

# Nouvelle frontière efficiente : la Capital Allocation Line !
plot(Z2[, "sd"], Z2[, "mean"], type = "b")

# ========================================================================== #
# LONG SHORT !
# ========================================================================== #

R <- colMeans(X)
S <- cov(X)
A <- cbind(1, diag(1, 30, 30), diag(-1, 30, 30))
B <- matrix(c(0, rep(-.2, 30), rep(-.2, 30)))

res <- solve.QP(S, .01*R, A, B, 1)

w <- matrix(res$solution, 30, 1)
rownames(w) <- Sym

po <- X %*% w
sd(po) * sqrt(260)

Chomolungma

Voici un exemple de chose que l’on peut faire avec quelques données topographiques et la fonction persp.

rm(list=ls())
setInternet2(TRUE)
url <- "collez ce lien entre les guillemets"
source(url)
# Charge *height*, une matrice [150, 150]
x <- 1:nrow(height)
y <- 1:ncol(height)
persp(x, y, height, zlim = c(0, 9e3), col = "lightsteelblue2", border = NA,
      shade = .2, phi = 25, theta = 30, axes = FALSE)

Et voilà :

Au cas où vous vous poseriez la question : oui, ce sont de vraies données ; ce que vous voyez ci-dessus, c’est le mont Everest en 3D.

Notez que comme mon zlim = c(0, 9e3) je suggère, la hauteur du graphique va du niveau de la mer à 9000 mètres.

Et si vous voulez vraiment vous faire plaisir (tout en découvrant à quoi servent theta et phi), essayez ça :

theta <- seq(30, 390, 5)
phi <- seq(0, 30, len = length(theta))

for(i in 1:length(theta)) {
      persp(x, y, height, zlim = c(0, 9e3), col = "lightsteelblue2",
            border = NA, shade = .2, phi = phi[i], theta = theta[i],
            axes = FALSE)
}

--- Edit

Après avoir installé le package rgl (install.packages("rgl")), exécutez ce qui suit (on suppose de height, x et y sont en mémoire) :

mx <- max(height)
xx <- drop(which(height == mx, TRUE, FALSE))
persp3d(x, y, height, col= "lightsteelblue2", zlim = c(0, 9e3), axes = FALSE,
      xlab = "", ylab = "", zlab = "")
points3d(xx[1], xx[2], mx, col = "red")
text3d(xx[1], xx[2], mx, "8840m", adj = c(.5, -1))

Vous obtiendrez le même type de graphique... mais interactif cette fois.

L'Everest, face du Kangshung, arrête Nord-Est, face Nord.

RBloomberg

Il semble que le package RBloomberg ne soit plus mis à jour. Je dis bien « il semble » parce qu’à titre personnel, je ne l’ai jamais utilisé.

Pour ceux que ça intéresse, voici mes fonctions : c’est un poil limité (pas de données intraday), on peut coder ça de manière beaucoup plus élégante mais ces petites fonctions offrent l'avantage de leur remarquable robustesse.

# ========================================================================= #
# BLOO.R
# ========================================================================= #

# ------------------------------------------------------------------------- #
# Options & libraries
if(! "RDCOMClient" %in% .packages()) library("RDCOMClient")

# ------------------------------------------------------------------------- #
# Utilities

# Converts R Date objets into COMDate objets.
b.BloombergDate = function(x) {
      res <- as.numeric(as.Date(x)) + 25569
      class(res) <- "COMDate"
      return(res)
}

# Converts COMDate objets into R Date objets.
b.RDate = function(x) {
      res <- as.Date("1899-12-30") + as.numeric(x)
      return(res)
}

# Replaces #N/A by NA.
b.NA = function(x) {
      x[substr(x, 1, 4) == "#N/A"] <- NA
      return(x)
}

# Turns a vector x into a list of vectors of max length n.
b.Frag = function(x, n) {
      nx <- length(x)
      v <- 1 + c(0, n * (1:ceiling(nx/n)))
      res <- lapply(1:(length(v)-1),
            function(i, x, v, nx) { x[v[i]:min(c(v[i+1]-1), nx)] },
            x = x, v = v, nx = nx)
      return(res)
}

# Error manager
is.error = function(x) class(x) == "try-error"

# ------------------------------------------------------------------------- #
# b.History
# Retrieves times series and forces them into a matrix.

# Args:
# tk         a vector of Bloomberg tickers
# fld        1 Bloomberg field, e.g. *LAST_PRICE* or *VOLUME* (etc...)
# start      start date (ISO 8601 format)
# end        end date (ISO 8601 format), defaults to Sys.Date()
# suffix     1 suffix (*Equity*, *Index* etc...)

# Output:
# A matrix with dates as row names and tk (without suffix) as colnames.

b.History = function(tk, fld, start, end = Sys.Date(), suffix = "") {
      if(length(fld) != 1) stop("'fld' must be of length 1")
      
      v <- seq.Date(as.Date(start), as.Date(end), by = 1)
      res <- matrix(NA, length(v), length(tk),, list(as.character(v), tk))
      
      b <- try(
            COMCreate("Bloomberg.Data.1"),
            silent = TRUE)
      
      if(! is.error(b)) {
            
            b[["Timeout"]] <- 12000
            b[["DisplayNonTradingDays"]] <- 64
            b[["NonTradingDayValue"]] <- 256
            b[["Periodicity"]] <- 1
            class(b) <- "COMIDispatch"
            
            d0 <- b.BloombergDate(start)
            d1 <- b.BloombergDate(end)
            
            U <- b.Frag(1:length(tk), 10)
            
            for(i in 1:length(U)) {
                  
                  tki <- paste(tk[U[[i]]], suffix)
                  
                  dta <- try(
                        b$BLPGetHistoricalData(
                              Security = tki,
                              Fields = fld,
                              StartDate = d0,
                              EndDate = d1),
                        silent = TRUE)
                  
                  if(! is.error(dta)) {
                        
                        d <- try(
                              lapply(dta[[1]],
                                    function(x) b.RDate(b.NA(unlist(x)))),
                              silent = TRUE)
                        r <- try(
                              lapply(dta[[2]],
                                    function(x) as.numeric(b.NA(unlist(x)))),
                              silent = TRUE)
                                          
                        if(! is.error(d) & ! is.error(r)) {
                              
                              tmp <- try(
                                    all(d[[1]] == do.call(cbind, d)),
                                    silent = TRUE)
                              chk <- if(! is.error(tmp)) tmp else FALSE
                              
                              if(chk) {
                                    ir <- as.character(d[[1]])
                                    ic <- tk[U[[i]]]
                                    res[ir, ic] <- do.call(cbind, r)
                              } else {
                                    for(j in 1:length(d)) {
                                          ir <- as.character(d[[j]])
                                          ic <- tk[U[[i]]][j]
                                          res[ir, ic] <- r[[j]]
                                    }
                              }
                        }
                  }
            }
      }
            
      rm(b)
      gc(verbose = FALSE)
      
      ix <- ! weekdays(v) %in% weekdays(as.Date("2000-01-01") + 0:1)
      res <- res[ix, ]
      return(res)
}

# ------------------------------------------------------------------------- #
# b.Point
# Retrieves one data for as many tickers you want.

# Args:
# tk         a vector of Bloomberg tickers
# f          1 Bloomberg field, e.g. *NAME* or *ID_ISIN* (etc...)
# string     logical: should output be coerced as string?

# Output:
# A vector of length(tk).

b.Point = function(tk, f, string = TRUE) {
      
      res <- rep(NA, length(tk))
      res <- if(string) as.character(res) else as.numeric(res)
      
      b <- try(
            COMCreate("Bloomberg.Data.1"),
            silent = TRUE)
      
      if(! is.error(b)) {
            
            b[["Timeout"]] <- 12000
            b[["DisplayNonTradingDays"]] <- 64
            b[["NonTradingDayValue"]] <- 256
            b[["Periodicity"]] <- 1
            class(b) <- "COMIDispatch"
            
            U <- b.Frag(1:length(tk), 300)
            
            for(i in 1:length(U)) {
                  
                  tki <- paste(tk[U[[i]]], "Equity")
                  
                  dta <- try(
                        b$BlpSubscribe(tki, f),
                        silent = TRUE)
                  
                  if(! is.error(dta)) {
                        tmp <- b.NA(unlist(dta))
                        if(string) {
                              res[U[[i]]] <- as.character(tmp)
                        } else {
                              res[U[[i]]] <- as.numeric(tmp)
                        }
                  }
            }
      }
      
      rm(b)
      gc(verbose = FALSE)
      
      return(res)
}

# ------------------------------------------------------------------------- #
# b.List
# Equivalent of b.Point but designed to retrieve index members. 

# Args:
# tk         1 Bloomberg tickers
# fld        1 Bloomberg field, defaults to *INDX_MEMBERS*
# suffix     1 suffix, defaults to *Index*

# Output:
# A vector.

b.List = function(tk, f = "INDX_MEMBERS", suffix = "Index") {
      
      res <- NULL
      
      b <- try(COMCreate("Bloomberg.Data.1"), silent = TRUE)
      
      if(! is.error(b)) {
            b[["Timeout"]] <- 12000
            b[["DisplayNonTradingDays"]] <- 64
            b[["NonTradingDayValue"]] <- 256
            b[["Periodicity"]] <- 1
            
            dta <- try(
                  b$BlpSubscribe(
                        Security = paste(tk, suffix),
                        Fields = f),
                  silent = TRUE)
            
            if(! is.error(dta)) {
                  res <- as.character(unlist(dta))
            } 
      }
      
      rm(b)
      gc(verbose = FALSE)
      
      return(res)
}
# -------------------------------- The End -------------------------------- #

Par exemple, pour récupérer la liste des composants du S&P 500 :

tk <- b.List("SPX")

Pour récupérer les noms :

b.Point(tk, "NAME", string = TRUE)

Pour les cours de clôture :

b.History(tk, "LAST_PRICE", from, to, suffix = "Index")

NB : il n'est pas impossible que j'ai laissé trainer quelques fonctions personnelles là-dedans (j'ai failli oublier de vous donner is.error). Vos retours d’expérience sont les bienvenus !

yHist

Une fonction pour télécharger des séries de prix via l’API de finance.yahoo.com.

Arguments :

  • sym, le symbole identifiant la série, i.e. YHOO pour les cours de l’actions de Yahoo! ;
  • from, date de début : un objet Date ou une chaine de caractère au format ISO 8601 (yyyy-mm-dd) ;
  • to, date de fin : un objet Date ou une chaine de caractère qui respecte la norme ISO 8601 (Sys.Date() par défaut) ;
  • freq, la fréquence des données : d pour des donnés quotidiennes, w pour des données hebdomadaires ou m pour des données mensuelles ;
  • flds, les champs à retourner en plus de Sym et Date symbolisés par une lettre : o pour l’ouverture, h pour le plus haut, l pour le plus bas, c pour la clôture, v pour le volume et a pour la clôture ajustée (tous par défaut).

yHist renvoie un dataframe. Utilisée avec plusieurs symboles, les données sont empilées les unes au-dessus des autres.

yHist = function(sym, from, to = Sys.Date(), freq = "d", flds = "ohlcva") {
      f <- as.Date(from[1])
      t <- as.Date(to[1])
      q <- tolower(substr(freq[1], 1, 1))
      if(! q %in% c("d", "w", "m")) stop("! freq %in% c('d', 'w', 'm')")
      tmp <- strsplit(flds[1], "")[[1]]
      o <- c(o = "Open", h = "High", l = "Low", c = "Close", v = "Volume",
            a="Adj.Close")
      if(any(! tmp %in% names(o))) stop("flds must be in 'ohlcva'")
      cn <- c("Sym", "Date", unname(o[tmp]))
      if(length(sym) > 1) {
            k <- lapply(sym, yHist, from = from, to = to, freq = freq,
                 flds = flds)
            res <- do.call(rbind, k)
      } else {
            url <- paste(c(
                 "http://ichart.yahoo.com/table.csv?",
                  "s=", sym,
                  "&a=", as.integer(format(f, "%m")) - 1,
                  "&b=", as.integer(format(f, "%d")),
                  "&c=", as.integer(format(f, "%Y")),
                  "&d=", as.integer(format(t, "%m")) - 1,
                  "&e=", as.integer(format(t, "%d")),
                  "&f=", as.integer(format(t, "%Y")),
                  "&g=", q,
                  "&ignore=.csv"), collapse = "")
            ans <- try(read.csv(url), silent = TRUE)
            if(class(ans) == "try-error") {
                  res <- NULL
            } else {
                  ans <- ans[order(ans$Date), ]
                  res <- subset(cbind(Sym = sym, ans), select = cn)
                  rownames(res) <- NULL
            }
      }
      return(res)
}

Exemple :

> yHist(c("GOOG", "AAPL", "YHOO", "MSFT"), from = "2014-01-01")
     Sym       Date    Open    High     Low   Close   Volume Adj.Close
1   GOOG 2014-01-02 1115.46 1117.75 1108.26 1113.12  1821400   1113.12
2   GOOG 2014-01-03 1115.00 1116.93 1104.93 1105.00  1666700   1105.00
3   GOOG 2014-01-06 1113.01 1118.86 1106.44 1117.32  1769300   1117.32
4   GOOG 2014-01-07 1125.00 1139.69 1121.16 1138.86  2552600   1138.86
5   GOOG 2014-01-08 1146.00 1147.32 1133.29 1141.23  2242500   1141.23
6   GOOG 2014-01-09 1143.44 1144.22 1125.56 1130.24  2084500   1130.24
...

Challenge #3

La Federal Reserve — la banque centrale des États-Unis — vous propose de télécharger un fichier .csv contenant les taux historiques des obligations d’État américaines.

Votre mission, si vous l’acceptez, consiste à écrire un script visant à stocker ces données dans un dataframe puis, puisque nous y sommes, à stocker ledit dataframe dans un fichier .RData que vous collerez où vous voulez.

La solution est ici.

Organisation

(NB : l'auteur s'est quelque peu laisser emporter sur ce post... à refaire...)

Coder c’est s’organiser. Plutôt que de me lancer dans explications sans fin, je vais vous donner un exemple d’organisation et vous laisser juger par vous-même de l’intérêt de la chose.

Voilà : pour des raisons qui me sont propres, je souhaite suivre les performances d’un certain nombre de fonds d’investissement. En l’occurrence, ce sont des mutual funds américains investis en actions étasuniennes que je compare — très logiquement — à l’indice S&P 500.

Ce que je voudrais, c’est disposer d’un historique complet de ces fonds (et du S&P 500) depuis le 31 mai 2012 (parce que !) ; historique prêt à l’emploi quand le besoin m’en prend.

Nous allons commencer, si vous le voulez bien, par créer un répertoire dédié à ce projet : mettez ça où vous voulez et appelez-le funds. Nous allons considérer que ce répertoire est situé à l’adresse :

C:/…/funds

Il se trouve que Yahoo! nous propose les données dont j’ai besoin à titre gracieux et que, mieux encore, l’entreprise dirigée par la très jolie Marissa Mayer a eu l’excellente idée de créer un système de requêtes en ligne qui génère des fichiers .csv à la demande.

Par exemple, si vous cliquez sur ce lien, vous récupèrerez les cours de l’indice S&P 500 du 31 mai 2012 au 31 mai 2013 au format .csv mais si vous cliquez sur celui-là, vous obtiendrez les cours d’Apple du 31 mai 2011 au 31 mai 2013.

En inspectant la manière dont ces liens sont construit (et, au besoin, en vous référant à cette page), vous pouvez assez facilement créer une fonction sous R qui permet de récupérer ces données pour n’importe quel titre et entre n’importe quelles dates.

Soit la fonction yget :

yget <- function(symb, from, to = Sys.Date(), quiet = TRUE) {
      
      # Construction de l'URL :
      url <- paste("http://ichart.finance.yahoo.com/table.csv?s=", 
            symb, 
            "&a=", as.numeric(format(from, "%m")) - 1, 
            "&b=", format(from, "%d"), 
            "&c=", format(from, "%Y"), 
            "&d=", as.numeric(format(to, "%m")) - 1, 
            "&e=", format(to, "%d"), 
            "&f=", format(to, "%Y"), 
            "&g=d&ignore=.csv", sep = "")
      
      # Lecture des données :
      res <- try(read.csv(url), silent = TRUE)
      
      # Notez l'usage de la fonction *try* qui renvoie un objet de classe
      # *try-error* en cas de problème :
      if(class(res) == "try-error") {
            
            # En cas d'erreur on renvoie un dataframe vide :
            res <- data.frame(Sym = NULL, Date = NULL, Open = NULL,
                  High = NULL, Low = NULL, Close = NULL, Adj.Close = NULL)
            
            # Et si *quiet = FALSE* on rajoute un avertissement :
            if(! quiet) warning("No data for ", sQuote(symb))
      
      # Sinon (i.e. si tout va bien) :
      } else {
            
            # On réordonne les données...
            res <- cbind(Sym = symb, res[order(res$Date), ])
      }
      
      # ... et on renvoie le résultat :
      return(res)
}

Ce qui nous donne :

> ans <- yget("^GSPC", as.Date("2012-05-31"))
> head(ans, 10)
      Sym       Date    Open    High     Low   Close     Volume Adj.Close
456 ^GSPC 2012-05-31 1313.09 1319.74 1298.90 1310.33 4557620000   1310.33
455 ^GSPC 2012-06-01 1309.87 1309.87 1277.25 1278.04 4669350000   1278.04
454 ^GSPC 2012-06-04 1278.29 1282.55 1266.74 1278.18 4011960000   1278.18
453 ^GSPC 2012-06-05 1277.82 1287.62 1274.16 1285.50 3403230000   1285.50
452 ^GSPC 2012-06-06 1285.61 1315.13 1285.61 1315.13 4268360000   1315.13
451 ^GSPC 2012-06-07 1316.15 1329.05 1312.68 1314.99 4258140000   1314.99
450 ^GSPC 2012-06-08 1314.99 1325.81 1307.77 1325.66 3497190000   1325.66
449 ^GSPC 2012-06-11 1325.72 1335.52 1307.73 1308.93 3537530000   1308.93
448 ^GSPC 2012-06-12 1309.40 1324.31 1306.62 1324.18 3442920000   1324.18
447 ^GSPC 2012-06-13 1324.02 1327.28 1310.51 1314.88 3506510000   1314.88

Formidable ! Copiez cette fonction et collez la dans un nouveau script que vous allez nommer ylib.R est sauvegarder à l’adresse :

C:/…/funds/ylib.R

L’air de rien, vous venez de poser le premier jalon de ce qui constituera — si ce genre de données vous intéressent — votre future de bibliothèque de fonctions dédiées à la récupération de données sur finance.yahoo.com.

Fermez ylib.R et ouvrez un nouveau script dans lequel nous allons écrire le code nécessaire à la mise à jour de nos données :

# ===========================================================================
# Yahoo! 
# ===========================================================================

# Le vicomte : Qu'est-ce que c'est que ça, s'il vous plaît ?
# Cyrano : C'est le titre.

# Ménage (au cas où) :
rm(list = ls(all = TRUE))

# Pointeur vers votre répertoire (vous n'avez que ça à taper !) :
.mydir <- "C:/.../funds"

# On charge la bibliothèque de fonctions :
source(file.path(.mydir, "ylib.R"))

# ---------------------------------------------------------------------------
# Params
# ---------------------------------------------------------------------------

# Le vicomte : Qu'est-ce que c'est que ça, s'il vous plaît ?
# Cyrano : Les paramètres de mise à jour.

# Nous voulons des données à partir du 31 mai 2012 :
from <- as.Date("2012-05-31")

# Une liste avec le S&P 500 et les fonds qui nous intéressent :
PeerGroup <- list(
    "^GSPC" = "S&P 500 Index",                           # (*)
      LMOPX = "Legg Mason Opportunity Trust",            # Bill Miller
      CPOAX = "Morgan Stanley Multi Cap Growth A",       # Lynch & team
      VCVLX = "Vanguard Capital Value Inv",              # Higgins & Palmer
      POAGX = "PrimeCap Odyssey Aggressive Growth Fund", # PrimeCap Team
      MXXVX = "Matthew 25"                               # Mark Mulholland
)

# (*) Ces types sont très très forts !

# ---------------------------------------------------------------------------
# Script 
# ---------------------------------------------------------------------------

# Le vicomte : Qu'est-ce que c'est que ça, s'il vous plaît ?
# Cyrano : Comment vous dire (sans être grossier) ?

# Téléchargement :
X <- lapply(names(PeerGroup), yget, from = from)

# Je construis un gros dataframe avec ce qui nous intéresse * :
Dta <- do.call(rbind, X)[, c("Sym", "Date", "Close", "Adj.Close")]
colnames(Dta) <- c("Sym", "Date", "NAV", "TR")
rownames(Dta) <- NULL

# (*) *Close* c'est le prix de clôture ; dans le cas d'un fonds on parlera
# plutôt de la *Net Asset Value* (NAV) et nous remplaçons *Adj.Close* (prix
# ajusté des dividendes) par *TR* (pour *Total Return*).

# On sauvegarde :
dump("Dta", file.path(.mydir, "Dta.R"))

Voilà un script bien ficelé dans lequel nous avons pris le soin de rajouter des commentaires (devant les #) afin de le rendre lisible à d’autres utilisateurs ou, d’ailleurs, à nous-mêmes si nous étions amenés à y revenir plus tard.

Enregistrons-le dans notre répertoire :

C:/…/funds/Update.R

Si vous exécutez ce script, vous constaterez qu’il créé un fichier Dta.R dans votre répertoire ; fichier dans lequel vous trouverez notre dataframe Dta. Il fera ça à chaque fois que vous l’exécuterez.

Seulement voilà : à ce stade, vous êtes encore obligés de lancer Rgui, d’ouvrir Update.R et de l’exécuter. C’est un peu fatiguant ; surtout si vous devez faire ça tous les jours. Nous allons donc nous simplifier le travail :

Avec un éditeur de texte — façon Bloc-notes de Windows — ouvrez un nouveau document texte et saisissez-y ceci :

"C:\Program Files\R\R-X.X.X\bin\Rscript.exe" "C:\...\funds\Update.R"

(R-X.X.X correspondant à votre version de R et "C:\...\funds\Update.R" étant l'adresse de Update.R.)

Lorsque c’est fait, sauvegardez ce fichier dans votre répertoire (appelez-le Update par exemple) en prenant soin de changer le suffixe .txt en .cmd :

Ce qui a pour effet de créer une chose de ce genre dans votre répertoire :

Double-cliquez sur la chose et vous obtenez :

... qui signifie que R est en train d'exécuter votre script.

Si c'est encore trop, rendez-vous dans le panneau de configuration, chapitre outils d’administration, où vous trouverez une section planificateur de tâches qui vous permettra de programmer le lancement automatiquement Update.cmd quand vous le souhaitez.

Voilà, ça n’est qu’un exemple et vous l’adapterez comme vous le souhaitez mais il a au moins l’avantage de vous montrer comment on se débarrasse d’une tâche quotidienne rébarbative au possible avec un peu d’organisation et quelques lignes de code.

Ce qui nous donne l’occasion de nous entraîner un peu…

La fonction f

Soit une matrice :

M <- matrix(head(rivers, 140), 28, 5)

Une chose que nous pourrions vouloir faire c'est calculer pour chaque ligne de M l'exposant a de la somme des éléments supérieurs à n.

En principe, on peut écrire ça de cette manière :

apply(X, 1, function(x, n, a) sum(x[x > n])^a, n = 350, a = 2)

Mais vous faites peut-être partie de ceux qui trouvent la syntaxe function(x, n, a)... est un peu verbeuse. Si c'est la cas, vous pouvez utiliser la fonction f :

f <- function(exp) {
      w <- all.names(substitute(exp), FALSE)
      if(any(w == "x")) w <- c("x", w[w != "x"])
      fo <- alist()
      length(fo) <- length(w) + 1
      names(fo) <- c(w, "call")
      fo[1:length(w)] <- expression(quote())
      fo$call <- substitute(exp)
      res <- as.function(fo)
      res
}

Cet animal-là est conçu pour transformer une expression en fonction anonyme avec pour particularité que si l'une des variables de l'expression exp s'appelle x, alors elle est nécessairement le premier argument de la fonction renvoyées par f.

Par exemple :

> foo <- f((x-y)^2)
> foo(3, 1)
[1] 4

Ce qui permet de raccourcir sensiblement notre appel à apply :

apply(X, 1, f(sum(x[x > n])^a), n = 350, a = 2)

Ou de faire des choses de ce genre :

> f(x+y)(1, 2)
[1] 3

Amusant non ?

Environnements #2

Lors d'un épisode précédent, nous avons compris comment fonctionnent les tiroirs de notre établi. Cette fois-ci, nous allons créer des boîtes au-dessus de .GlobalEnv.

Soit deux vecteurs...

rm(list = ls(all = TRUE))
y <- x <- 1

... et une boîte :

box <- new.env()

Sans surprise, on constate que :

> ls()
[1] "box" "x"   "y"

Par ailleurs :

> class(box)
[1] "environment"
> ls(box)
character(0)
> parent.env(box)

C'est-à-dire que notre boîte box est un environnement vide et parent.env nous précise qu'elle est posée juste au-dessus de .GlobalEnv (contrairement aux tiroirs qui, souvenez-vous, sont en-dessous).

Nous allons maintenant ranger deux objets x et y dans notre boîte :

box$x <- 1:3
box$y <- 2

À la surface de notre établi, on a toujours :

> ls()
[1] "box" "x"   "y"

Mais dans la boîte :

> ls(box)
[1] "x" "y"

Avec :

> get("x", box)
[1] 1 2 3

Et (autre méthode) :

> box$y
[1] 2

Si nous sommons le x et le y posés sur .GlobalEnv :

> x+y
[1] 4

Mais si nous évaluons l'expression x + y dans box :

> eval(expression(x+y), envir = box)
[1] 3 4 5

Logique n'est-ce pas ?

Bien. Maintenant, considérez cette fonction :

foo <- function(x, y) {
      res <- x + y^2
      return(res)
}

Si nous l'exécutons dans .GlobalEnv, nous obtenons :

> foo(x, y)
[1] 6

Question : pourquoi res n’apparaît-il pas :

> ls()
[1] "box" "foo" "x"   "y"

Eh bien c'est fort simple : parce que foo a évalué res dans un environnement distinct de .GlobalEnv. On peut le vérifier comme suit :

foo <- function(x, y) {
      res <- x + y^2
      e <- environment()
      return(e)
}

Ce qui nous donne :

> ans <- foo(x, y)
> ans
<environment: 0x105d9c84>
> ls(ans)
[1] "e"   "res" "x"   "y"
> ans$res
[1] 6

Par ailleurs...

> parent.env(ans)
<environment: R_GlobalEnv>

... nous signale que cet environnement est posé au-dessus de .GlobalEnv.

Ça n’est pas du tout trivial. Voici pourquoi :

foo <- function(x) x + a

On peut vérifier que :

> a <- 5
> foo(x)
[1] 7

Tandis qu'avec :

rm(list = "a")
box$a <- 5

Gruik !

> foo(x)
Erreur dans foo(x) : objet 'a' introuvable

R a cherché a dans l'environnement de foo, dans .GlobalEnv, dans tous les tiroirs jusqu'au tiroir vide et rien, rien de rien : pas la moindre trace de a ; d'où le message d'erreur. En revanche, R n'est pas allée voir dans box parce que box est un environnement parallèle à celui de foo.

Shématiquement, en notant <foo> l'environnement de foo, on peut représenter le chemin de recherche de R comme suit :

<foo> --> |
          | <.GlobalEnv> --> <rgl> --> (...) --> <base> --> <vide>
<box>     |
C’est pour cette raison que ce qui est évalué dans foo n’est pas accessible dans .GlobalEnv et que ce qui est dans box n’est pas accessible à foo.

En revanche, avec :

rm(list = ls(all = TRUE))
box <- new.env()
box$a <- 1
foo <- function(x) x + box$a

On obtient :

> foo(1)
[1] 2

Parce que box étant posé que .GlobalEnv, il est bien sur le chemin de recherche et il suffit de dire à R de remonter chercher a dans cet environnement pour que ça fonctionne.

Si vous souhaitez pousser des objets hors de la boîte d'une fonction, vous pouvez faire des choses comme ça :

rm(list = ls(all = TRUE))
check <- function() {
      if(exists(".count", envir = .GlobalEnv)) {
            .count <<- .count + 1
      } else {
            .count <<- 1
      }
      return(.count)
}

Ce qui donne :

> check()
[1] 1
> check()
[1] 2
> check()
[1] 3
> .count
[1] 3

Ou des choses comme ça :

rm(list = ls(all = TRUE))
make.rot <- function(n) {
 foo <- function(x) {
            ix <- (match(strsplit(x, "")[[1]], letters) + n) %% 26
            paste(letters[ix], collapse = "")
        }
 fname <- paste("rot", n, sep = "")
 assign(fname, foo, .GlobalEnv)
 if(exists(fname, .GlobalEnv)) TRUE else FALSE
}

D'où :

> make.rot(13)
[1] TRUE
> ls()
[1] "make.rot" "rot13"
> rot13("obawbhe")
[1] "bonjour"

Environnements #1

Un bon artisan range son établi et n’y laisse jamais — Ô grand jamais — trainer ses outils. C’est à ça que servent les tiroirs en dessous : à ranger des outils et même, éventuellement, un peu de matière première.

R fonctionne comme ça.

La surface de votre établi s’appelle .GlobalEnv ou, plus communément, espace de travail.

Pour inspecter la surface de .GlobalEnv et lister les objets qui se trouvent dessus (outils ou autres) vous utilisez la fonction ls sans argument ce qui revient, par défaut, à établir une liste des objets présents dans l’environnement dans lequel vous travaillez ; en l’occurrence, la surface de votre établi.

Lorsque vous démarrer une nouvelle cession — et à supposer que vous n’avez pas sauvegardé la cession précédents (ce qui signifie que vous avez sciemment laissé des trucs trainer sur votre établi depuis) — vous devriez obtenir :

> ls()
character(0)

Ce qui revient à dire :

> ls(.GlobalEnv)
character(0)

Bref, la surface de votre établi est absolument vide.

Commençons par y déposer un peu de matière première ; en l'espèce x et y, deux vecteurs numérique de longueur 1.

> y <- x <- 1
> x
[1] 1
> y
[1] 1

De telle sorte que :

> ls()
[1] "x" "y"

Sur ces deux objets, nous allons appeler la fonction c qui est notre tube de colle universelle :

> c(x, y)
[1] 1 1

À ce stade, vous pourriez vous demander d’où sort c puisque que ls ne l’a pas listée parmi les objets présents à la surface de notre établi. C’est une excellente question que vous auriez d’ailleurs pu vous poser pour ls aussi.

La réponse est très simple : vous avez trouvé ces deux outils dans un des tiroirs de votre établi. En l’occurrence l’avant dernier, celui sur lequel il est écrit : .BaseNamespaceEnv.

Voilà comment ça fonctionne : à chaque fois que vous avez besoin de votre tube de colle (la fonction c), vous commencez par chercher quelque chose qui porte l’étiquette c (et qui soit un outil, on y reviendra) à la surface de votre établi (.GlobalEnv). Si vous trouvez un outil étiqueté c sur .GlobalEnv, vous utilisez celui-là mais dans le cas contraire, vous cherchez dans le premier de vos tiroirs puis, dans le second, le troisième… jusqu’à ce que vous trouviez un outil (un objet de classe function) étiqueté c.

Normalement, la fonction c se trouve dans l’avant dernier tiroir — celui où il est écrit .BaseNamespaceEnv et qui contient tous vos outils de base — sachant que le dernier tiroir et toujours vide.

Illustrons ça avec un exemple :

Soit une nouvelle fonction c définie comme :

c <- function(...) {
      a <- list(...)
      sum(unlist(a) * 10^((length(a)-1):0))
}

Cette fonction est définie sur .GlobalEnv de telle sorte que quand vous appellerez c, c'est elle que vous utiliserez par défaut :

> c(1, 1)
[1] 11

Mais si vous souhaitez ponctuellement utiliser la fonction c qui se trouve dans .BaseNamespaceEnv (l'avant-dernier tiroir), vous pouvez le faire en écrivant :

> base::c(1, 1)
[1] 1 1

Dans R, la surface de votre établi et chacun de ses tiroirs est un objet de classe environment :

> class(.GlobalEnv)
[1] "environment"
> class(.BaseNamespaceEnv)
[1] "environment"

Si vous souhaitez lister les objets présents dans .BaseNamespaceEnv (toutes les fonctions du package base de R) :

ls(.BaseNamespaceEnv)

Ou :

ls(baseenv())

(C'est un gros tiroir...)

Si vous souhaitez savoir dans quel ordre R cherche dans ses tiroirs :

> search()
 [1] ".GlobalEnv"        "package:rgl"       "package:stats"    
 [4] "package:graphics"  "package:grDevices" "package:utils"    
 [7] "package:datasets"  "package:methods"   "Autoloads"        
[10] "package:base"

Comme spécifié ci-dessus, R va d'abord chercher c dans .GlobalEnv, puis dans package:rgl, puis dans package:stats ... jusqu'à package:base mais pas dans le dernier tiroir (emptyenv()) puisqu'il est toujours vide. En l'occurrence, comme vous avez créé une fonction c dans .GlobalEnv, il utilise celle-là.

Mais si vous la détruisez :

> rm(list = "c")
> c(1, 1)
[1] 1 1

C'est la fonction c de package:base qui reprend le dessus.

Lorsque vous exécutez :

> diff
function (x, ...) 
UseMethod("diff")
<bytecode: 0x07279760>
<environment: namespace:base>

La ligne <environment: namespace:base> vous signale que diff se trouve dans l'environnement de base (a.k.a. avant-dernier tiroir).

Deux petites précisions à ce stade :

Primo, R gère séparément les fonctions et les autres objets de telle sorte que vous pouvez, par exemple, appeler la fonction diff (base::diff) sur un objet diff créé dans votre espace de travail :

> diff <- seq(0, 10, 2)
> diff(diff)
[1] 2 2 2 2 2

Deuxio, ls() ne liste pas vraiment tous les objets en mémoire :

> rm(list=ls())
> x <- 1
> .coucou <- 2
> ls()
[1] "x"
> rm(list = ls())
> .coucou
[1] 2

Dans ce cas, utilisez plutôt :

> ls(all = TRUE)
[1] ".coucou"      ".Random.seed"

(Un petit >?.Random.seed répondra à vos questions.)

Bien. Maintenant que vous savez que votre espace de travail est un environnement, qu'il y a des tiroirs en dessous, que ces tiroirs s'appellent sont aussi des environnements et que R cherche les objets que vous utilisez dans vos tiroirs en suivant un ordre prédéfinis, nous allons créer des boîtes au-dessus de .GlobalEnv ; boîtes qui seront aussi des environnements ; suivez le guide...

Structure

Soit m, un objet :

m <- structure(
      1:9,
      .Dim = c(3L, 3L),
      .Dimnames = list(c("A", "B", "C"), c("Jan", "Feb", "Mar")))

Dans la console :

> m
  Jan Feb Mar
A   1   4   7
B   2   5   8
C   3   6   9
> class(m)
[1] "matrix"

C’est-à-dire que m est une matrice, exactement la même que sui vous aviez exécuté :

m <- matrix(1:9, 3, 3,, list(c("A", "B", "C"), c("Jan", "Feb", "Mar")))

Cette forme de m exprimée avec structure(...), c'est celle que vous obtiendriez en utilisant la fonction dump ; c'est une représentation de m avec ses données (1:9) et ses attributs (.Dim et .Dimnames).

> attributes(m)
$dim
[1] 3 3

$dimnames
$dimnames[[1]]
[1] "A" "B" "C"

$dimnames[[2]]
[1] "Jan" "Feb" "Mar"

Pour vous débarrasser de l'attribut .Dimnames, par exemple, il suffit d'écrire :

> attr(m, "dimnames") <- NULL
> m
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

La fonction unname fait la même chose. A ce stade, la représentation de votre objet m ressemble à ceci :

m <- structure(1:9, .Dim = c(3L, 3L))

De la même manière, vous pouvez casser les dimensions de m en exécutant :

> attr(m, "dim") <- NULL
> m
[1] 1 2 3 4 5 6 7 8 9

Si vous faites ça, m n'est plus que :

m <- 1:9

Et :

> class(m)
[1] "integer"

C'est exactement ce que vous faites lorsque vous exécutez une fonction comme as.numeric sur une matrice, vous cassez ses dimensions et, par la même occasion, supprimez les noms desdites dimensions.

Voici quelques structures-types :

# vector (with names)
v <- structure(
      c(100, 99, 101),
      .Names = c("Jan", "Feb", "Mar"))

# matrix (with dimnames)
m <- structure(
      1:9,
      .Dim = c(3L, 3L),
      .Dimnames = list(c("A", "B", "C"), c("Jan", "Feb", "Mar")))

# factor
f <- structure(
 c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L),
 .Label = c(month.abb),
 class = "factor")

# list
l <- structure(
      list( Month = c("Jan", "Feb", "Mar"),
            Value = c(100, 99, 101)),
      .Names = c("Month", "Value"))

# data.frame
d <- structure(
      list( Month = structure(
                  1:12,
             .Label = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul",
                  "Aug", "Sep", "Oct", "Nov", "Dec"),
             class = "factor"),
            Value = 1:12),
      .Names = c("Month", "Value"),
      row.names = c(NA, -3L),
      class = "data.frame")

Où vous comprenez pourquoi il ne sert à rien d'essayer de casser l'attribut .Dim d'un dataframe (il n'y a rien à casser). En revanche, vous pouvez faire ça :

> class(d) <- "list"
> d
$Month
 [1] Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
Levels: Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec

$Value
 [1]  1  2  3  4  5  6  7  8  9 10 11 12

attr(,"row.names")
[1] 1 2 3
> class(d)
[1] "list"

Méthodes

À toutes fins utiles je rappelle que, sur vecteur, la fonction rev fait ça :

> rev(1:3)
[1] 3 2 1

Maintenant, construisons une matrice :

rm(list=ls())
x <- matrix(1:9, 3, 3,, dimnames = list(LETTERS[1:3], month.abb[1:3]))

De telle sorte que :

> x
  Jan Feb Mar
A   1   4   7
B   2   5   8
C   3   6   9
> class(x)
[1] "matrix"
> mode(x)
[1] "numeric"

Si vous essayez rev sur x :

> rev(x)
[1] 9 8 7 6 5 4 3 2 1

Ce n'est pas ce que vous vouliez. Ce que vous auriez aimé, c'est que rev inverse l'ordre des lignes de votre matrice.

Bien sûr, vous pouvez créer une fonction ad hoc mais — si c’est vraiment ce que vous souhaitez — vous pouvez aussi modifier le comportement de rev sur les matrices. Voici comment faire :

rev.matrix <- function(x) x[nrow(x):1,, drop = FALSE]

Cette syntaxe generique.classe vous permet de créer une méthode pour rev sur la classe d'objet matrix. Une fois chargée en mémoire :

> rev(x)
  Jan Feb Mar
C   3   6   9
B   2   5   8
A   1   4   7

Et maintenant, créez une méthode pour rev sur la classe d'objet bidule qui fait la même chose mais sur l'ordre des colonnes :

rev.bidule <- function(x) x[, ncol(x):1, drop = FALSE]

Modifiez la classe de x :

class(x) <- "bidule"

Et testez :

> class(x)
[1] "bidule"
> mode(x)
[1] "numeric"
> rev(x)
  Mar Feb Jan
A   7   4   1
B   8   5   2
C   9   6   3

Fondamentalement, c'est à ça que servent les classes d'objets.

Last Observation Copied Forward

Si vous travaillez sur des séries temporelles, il y a une chose que vous avez toujours voulu faire : remplacer les données manquantes — i.e. les NA — par la dernière observation connue. C’est l’idée de la fonction locf pour Last Observation Copied Forward.

C’est-à-dire que qui vous avez :

> x <- c(1:3, NA, 5:6)
> x

[1] 1 2 3 NA 5 6

Vous voulez :

> locf(x)
[1] 1 2 3 3 5 6

La version la plus courte (et performante!) que j'ai trouvé s'écrit en une ligne :

x[c(NA, which(!is.na(x)))[cumsum( !is.na(x) ) + 1]]

Encapsulée dans une fonction récursive qui traite les matrice par colonnes, ça donne :

locf = function(x) {
        if( is.matrix(x) ) {
                res <- apply(x, 2, locf)
                dimnames(res) <- dimnames(x)
        } else {
                res <- x[c(NA, which(!is.na(x)))[cumsum( !is.na(x) ) + 1]]
                names(res) <- names(x)
        }
        return(res)
}

Palindromes

Avec :

url <- "collez ce lien entre les guillemets"
dico <- scan(url, what = "character", sep = "\n")

> length(dico)
[1] 336531

Exécutez ça :

sapply(dico, function(a) {
      paste(rev(strsplit(a, "")[[1]]), collapse = "")
} ) -> ocid

Ce qui nous donne une liste de palindromes:

> dico[which(dico == ocid)]
 [1] "a"         "à"         "alla"      "ana"       "ara"       "aviva"    
 [7] "axa"       "bob"       "cc"        "elle"      "erre"      "essayasse"
[13] "esse"      "été"       "étêté"     "eue"       "gag"       "ici"      
[19] "kayak"     "lebel"     "nanan"     "non"       "pep"       "pop"      
[25] "radar"     "ressasser" "retâter"   "rotor"     "sagas"     "salas"    
[31] "sanas"     "sapas"     "sas"       "sassas"    "selles"    "semâmes"  
[37] "sénés"     "sennes"    "serres"    "ses"       "sexes"     "shahs"    
[43] "sis"       "snobons"   "solos"     "sonos"     "sus"       "tâtât"    
[49] "tôt"       "tut"       "tût"       "y"

... ou une liste d'anacyclique :

> intersect(dico, ocid)
  [1] "a"         "à"         "adoré"     "ados"      "adulé"     "affins"   
  [7] "ah"        "ail"       "ailla"     "aimé"      "air"       "aira"     
 [13] "alevin"    "alla"      "allia"     "alliacé"   "amis"      "an"       
 [19] "ana"       "angor"     "animal"    "annotas"   "annotât"   "annoté"   
 [25] "ara"       "ares"      "aria"      "arum"      "as"        "assit"    
 [31] "aval"      "avaler"    "aviva"     "axa"       "bac"       "bob"      
 [37] "bol"       "bons"      "but"       "cab"       "cal"       "calf"     
 [43] "camus"     "cas"       "casser"    "cc"        "ces"       "cor"      
 [49] "coté"      "crut"      "dus"       "écaffe"     etc...