## Copyright (C) 1997-2003 Adrian Trapletti ## ## This program is free software; you can redistribute it and/or modify ## it under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 2, or (at your option) ## any later version. ## ## This program is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## A copy of the GNU General Public License is available via WWW at ## http://www.gnu.org/copyleft/gpl.html. You can also obtain it by ## writing to the Free Software Foundation, Inc., 59 Temple Place, ## Suite 330, Boston, MA 02111-1307 USA. ## ## Financial time series analysis ## portfolio.optim <- function (x, ...) UseMethod ("portfolio.optim") portfolio.optim.ts <- function (x, ...) { if(!is.ts(x)) stop("method is only for time series") if(NCOL(x) == 1) stop("x is not a multivariate time series") res <- portfolio.optim.default(as.matrix(x), ...) res$px <- ts(res$px, start = start(x), frequency = frequency(x)) return(res) } portfolio.optim.default <- function(x, pm = mean(x), riskless = FALSE, shorts = FALSE, rf = 0.0, reslow = NULL, reshigh = NULL, covmat = cov(x), ...) { if(NCOL(x) == 1) stop("x is not a matrix") if(any(is.na(x))) stop("NAs in x") k <- dim(x)[2] if(!is.matrix(covmat)) { stop("covmat is not a matrix") } if((dim(covmat)[1] !=k) || (dim(covmat)[2] !=k)) { stop("covmat has not the right dimension") } Dmat <- covmat dvec <- rep.int(0, k) big <- 1e+100 if(!is.null(reslow) && is.null(reshigh)) { reshigh <- rep.int(big, k) } if(is.null(reslow) && !is.null(reshigh)) { reslow <- -rep.int(big, k) } if(!is.null(reslow)) { if(!is.vector(reslow)) { stop("reslow is not a vector") } if(length(reslow) != k) { stop("reslow has not the right dimension") } } if(!is.null(reshigh)) { if(!is.vector(reshigh)) { stop("reshigh is not a vector") } if(length(reshigh) != k) { stop("reshigh has not the right dimension") } } if(riskless) { a1 <- colMeans(x) - rf if(shorts) { a2 <- NULL b2 <- NULL } else { a2 <- matrix(0, k, k) diag(a2) <- 1 b2 <- rep.int(0, k) } if(!is.null(reslow) && !is.null(reshigh)) { a3 <- matrix(0, k, k) diag(a3) <- 1 Amat <- t(rbind(a1, a2, a3, -a3)) b0 <- c(pm-rf, b2, reslow, -reshigh) } else { Amat <- t(rbind(a1, a2)) b0 <- c(pm-rf, b2) } res <- solve.QP(Dmat, dvec, Amat, bvec=b0, meq=1) } else { a1 <- rep.int(1, k) a2 <- colMeans(x) if(shorts) { if(!is.null(reslow) && !is.null(reshigh)) { a3 <- matrix(0, k, k) diag(a3) <- 1 Amat <- t(rbind(a1, a2, a3, -a3)) b0 <- c(1, pm, reslow, -reshigh) } else { Amat <- t(rbind(a1, a2)) b0 <- c(1, pm) } } else { a3 <- matrix(0, k, k) diag(a3) <- 1 b3 <- rep.int(0, k) if(!is.null(reslow) && !is.null(reshigh)) { Amat <- t(rbind(a1, a2, a3, a3, -a3)) b0 <- c(1, pm, b3, reslow, -reshigh) } else { Amat <- t(rbind(a1, a2, a3)) b0 <- c(1, pm, b3) } } res <- solve.QP(Dmat, dvec, Amat, bvec=b0, meq=2) } y <- c(tcrossprod(res$solution, x)) ans <- list(pw = res$solution, px = y, pm = mean(y), ps = sd(y)) return(ans) } get.hist.quote <- function (instrument = "^gdax", start, end, quote = c("Open", "High", "Low", "Close"), provider = c("yahoo", "oanda"), method = NULL, origin = "1899-12-30", compression = "d", retclass = c("zoo", "its", "ts"), quiet = FALSE, drop = FALSE) ## Added new argument 'compression'. ## May be "d", "w" or "m", for daily weekly or monthly. ## Defaults to "d". ## John Bollinger, 2004-10-27, www.BollingerBands.com, bbands@yahoo.com ## ## Changed POSIXct class to Date class, 2005-03-31 ## Added drop= arg, 2006-06-04 { if(missing(start)) start <- "1991-01-02" if(missing(end)) end <- format(Sys.Date() - 1, "%Y-%m-%d") provider <- match.arg(provider) retclass <- match.arg(retclass) start <- as.Date(start) end <- as.Date(end) if(is.null(method)) { method <- getOption("download.file.method") if(is.null(method)) method <- "auto" } if(provider == "yahoo") { url <- paste("http://chart.yahoo.com/table.csv?s=", instrument, format(start, paste("&a=", as.character(as.numeric(format(start, "%m"))-1), "&b=%d&c=%Y", sep = "")), format(end, paste("&d=", as.character(as.numeric(format(end, "%m"))-1), "&e=%d&f=%Y", sep = "")), "&g=", compression, "&q=q&y=0&z=", instrument, "&x=.csv", sep = "") destfile <- tempfile() ## Try downloading several times. ## Could make this controllable ... i <- 1L repeat { status <- tryCatch(download.file(url, destfile, method = method, quiet = quiet), error = identity) if(!inherits(status, "error") && (status == 0)) break unlink(destfile) if(i >= 5L) { if(inherits(status, "error")) stop(conditionMessage(status)) else stop(sprintf("download error, status %d", status)) } message("download error, retrying ...") i <- i + 1L } nlines <- length(count.fields(destfile, sep = "\n")) if(nlines == 1) { unlink(destfile) stop(paste("no data available for", instrument)) } ## Yahoo includes rows concerning dividends, ## hence need fill = TRUE and na.omit x <- read.table(destfile, header = TRUE, sep = ",", as.is = TRUE, fill = TRUE) x <- na.omit(x) ## As of 2012-07-07, ## get.hist.quote("ibm", start = "2012-07-06", end = "2012-07-06") ## gives two entries for 2012-07-06 ... ## Fix (hopefully only needed temporarily?) suggested by AZ. if(nrow(x) >= 2L && x[1L, 1L] == x[2L, 1L]) { warning("first date duplicated, first instance omitted") x <- x[-1L, , drop = FALSE] } ## Debug ## cat("read.table: start =", x[NROW(x),"Date"], "\n") ## cat("read.table: end =", x[1,"Date"], "\n") unlink(destfile) names(x) <- gsub("\\.", "", names(x)) nser <- pmatch(quote, names(x)[-1]) + 1 if(any(is.na(nser))) stop("this quote is not available") n <- nrow(x) ## Yahoo now seems to format dates as %Y-%m-%d. dat <- as.Date(as.character(x[, 1]), "%Y-%m-%d") if(!quiet && (dat[n] != start)) cat(format(dat[n], "time series starts %Y-%m-%d\n")) if(!quiet && (dat[1] != end)) cat(format(dat[1], "time series ends %Y-%m-%d\n")) if(retclass == "ts") { jdat <- unclass(julian(dat, origin = as.Date(origin))) ## We need unclass() because 1.7.0 does not allow adding a ## number to a "difftime" object. ind <- jdat - jdat[n] + 1 y <- matrix(NA, nrow = max(ind), ncol = length(nser)) y[ind, ] <- as.matrix(x[, nser, drop = FALSE]) colnames(y) <- names(x)[nser] y <- y[, seq_along(nser), drop = drop] return(ts(y, start = jdat[n], end = jdat[1])) } else { x <- as.matrix(x[, nser, drop = FALSE]) rownames(x) <- NULL y <- zoo(x, dat) y <- y[, seq_along(nser), drop = drop] if(retclass == "its") { if(inherits(tryCatch(getNamespace("its"), error = identity), "error")) warning("package its could not be loaded: zoo series returned") else { index(y) <- as.POSIXct(index(y)) y <- its::as.its(y) } } return(y) } } else if(provider == "oanda") { if(!missing(quote)) { warning("argument 'quote' ignored for provider 'oanda'") } if(!missing(compression)) { warning("argument 'compression' ignored for provider 'oanda'") } url <- paste("http://www.oanda.com/convert/fxhistory?lang=en&date1=", format(start, "%m"), "%2F", format(start, "%d"), "%2F", format(start, "%y"), "&date=", format(end, "%m"), "%2F", format(end, "%d"), "%2F", format(end, "%y"), "&date_fmt=us&exch=", unlist(strsplit(instrument, split = "/"))[1], "&exch2=&expr=", unlist(strsplit(instrument, split = "/"))[2], "&expr2=&margin_fixed=0&&SUBMIT=Get+Table&format=ASCII&redirected=1", sep="") destfile <- tempfile() status <- download.file(url, destfile, method = method, quiet = quiet) if(status != 0) { unlink(destfile) stop(paste("download error, status", status)) } x <- readLines(destfile, warn = quiet) unlink(destfile) if(length(grep("Sorry", x)) > 0) { msg <- unlist(strsplit(gsub("<[a-zA-Z0-9\\/]*>", "", x[grep("Sorry", x)]), split = " ")) msg <- paste(msg[msg != ""], collapse = " ") stop("Message from Oanda: ", msg) } ## This used to have ## first <- which(substr(x, 1, 5) == "
") ## last <- which(x == "") - 1 ## but at least on 2007-11-04 we had an instance of
not ## being at the beginning of a line andnot being on a ## line of its own ... first <- grep("
", x, fixed = TRUE) last <- grep("", x, fixed = TRUE) - 1 if((length(first) != 1) || (length(last) != 1) || (last < first)) { stop(paste("no data available for", instrument)) } x[first] <- sub(".*
", "", x[first]) con <- textConnection(x[first:last]) on.exit(close(con)) x <- scan(con, what = list(character(), double()), quiet = TRUE) dat <- as.Date(x[[1]], format = "%m/%d/%Y") n <- length(dat) if(!quiet && (dat[1] != start)) cat(format(dat[1], "time series starts %Y-%m-%d\n")) if(!quiet && (dat[n] != end)) cat(format(dat[n], "time series ends %Y-%m-%d\n")) if(retclass == "ts") { jdat <- unclass(julian(dat, origin = as.Date(origin))) ind <- jdat - jdat[1] + 1 y <- rep.int(NA, max(ind)) y[ind] <- x[[2]] return(ts(y, start = jdat[1], end = jdat[n])) } else { y <- zoo(x[[2]], dat) if(retclass == "its") { if(inherits(tryCatch(getNamespace("its"), error = identity), "error")) warning("package its could not be loaded: zoo series returned") else { index(y) <- as.POSIXct(index(y)) y <- its::as.its(y) } } return(y) } } else stop("provider not implemented") } maxdrawdown <- function(x) { if(NCOL(x) > 1) stop("x is not a vector or univariate time series") if(any(is.na(x))) stop("NAs in x") cmaxx <- cummax(x)-x mdd <- max(cmaxx) to <- which(mdd == cmaxx) from <- double(NROW(to)) for (i in 1:NROW(to)) from[i] <- max(which(cmaxx[1:to[i]] == 0)) return(list(maxdrawdown = mdd, from = from, to = to)) } sharpe <- function(x, r = 0, scale = sqrt(250)) { if(NCOL(x) > 1) stop("x is not a vector or univariate time series") if(any(is.na(x))) stop("NAs in x") if(NROW(x) == 1) return(NA) else { y <- diff(x) return(scale * (mean(y)-r)/sd(y)) } } sterling <- function(x) { if(NCOL(x) > 1) stop("x is not a vector or univariate time series") if(any(is.na(x))) stop("NAs in x") if(NROW(x) == 1) return(NA) else { return((x[NROW(x)]-x[1]) / maxdrawdown(x)$maxdrawdown) } } plotOHLC <- function(x, xlim = NULL, ylim = NULL, xlab = "Time", ylab, col = par("col"), bg = par("bg"), axes = TRUE, frame.plot = axes, ann = par("ann"), main = NULL, date = c("calendar", "julian"), format = "%Y-%m-%d", origin = "1899-12-30", ...) { if ((!is.mts(x)) || (colnames(x)[1] != "Open") || (colnames(x)[2] != "High") || (colnames(x)[3] != "Low") || (colnames(x)[4] != "Close")) stop("x is not a open/high/low/close time series") xlabel <- if (!missing(x)) deparse(substitute(x)) else NULL if (missing(ylab)) ylab <- xlabel date <- match.arg(date) time.x <- time(x) dt <- min(lag(time.x)-time.x)/3 if (is.null(xlim)) xlim <- range(time.x) if (is.null(ylim)) ylim <- range(x[is.finite(x)]) plot.new() plot.window(xlim, ylim, ...) segments(time.x, x[, "High"], time.x, x[, "Low"], col = col[1], bg = bg) segments(time.x - dt, x[, "Open"], time.x, x[, "Open"], col = col[1], bg = bg) segments(time.x, x[, "Close"], time.x + dt, x[, "Close"], col = col[1], bg = bg) if (ann) title(main = main, xlab = xlab, ylab = ylab, ...) if (axes) { if (date == "julian") { axis(1, ...) axis(2, ...) } else { n <- NROW(x) lab.ind <- round(seq(1, n, length=5)) D <- as.vector(time.x[lab.ind]*86400) + as.POSIXct(origin, tz = "GMT") DD <- format.POSIXct(D, format = format, tz ="GMT") axis(1, at=time.x[lab.ind], labels=DD, ...) axis(2, ...) } } if (frame.plot) box(...) }