1## Copyright (C) 1997-2003 Adrian Trapletti 2## 3## This program is free software; you can redistribute it and/or modify 4## it under the terms of the GNU General Public License as published by 5## the Free Software Foundation; either version 2, or (at your option) 6## any later version. 7## 8## This program is distributed in the hope that it will be useful, but 9## WITHOUT ANY WARRANTY; without even the implied warranty of 10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 11## General Public License for more details. 12## 13## A copy of the GNU General Public License is available via WWW at 14## http://www.gnu.org/copyleft/gpl.html. You can also obtain it by 15## writing to the Free Software Foundation, Inc., 59 Temple Place, 16## Suite 330, Boston, MA 02111-1307 USA. 17 18## 19## Financial time series analysis 20## 21 22portfolio.optim <- function (x, ...) UseMethod ("portfolio.optim") 23 24portfolio.optim.ts <- 25function (x, ...) 26{ 27 if(!is.ts(x)) 28 stop("method is only for time series") 29 if(NCOL(x) == 1) 30 stop("x is not a multivariate time series") 31 res <- portfolio.optim.default(as.matrix(x), ...) 32 res$px <- ts(res$px, start = start(x), frequency = frequency(x)) 33 return(res) 34} 35 36portfolio.optim.default <- 37function(x, pm = mean(x), riskless = FALSE, shorts = FALSE, 38 rf = 0.0, reslow = NULL, reshigh = NULL, covmat = cov(x), ...) 39{ 40 if(NCOL(x) == 1) 41 stop("x is not a matrix") 42 if(any(is.na(x))) 43 stop("NAs in x") 44 k <- dim(x)[2] 45 if(!is.matrix(covmat)) { 46 stop("covmat is not a matrix") 47 } 48 if((dim(covmat)[1] !=k) || (dim(covmat)[2] !=k)) { 49 stop("covmat has not the right dimension") 50 } 51 Dmat <- covmat 52 dvec <- rep.int(0, k) 53 big <- 1e+100 54 if(!is.null(reslow) && is.null(reshigh)) { 55 reshigh <- rep.int(big, k) 56 } 57 if(is.null(reslow) && !is.null(reshigh)) { 58 reslow <- -rep.int(big, k) 59 } 60 if(!is.null(reslow)) { 61 if(!is.vector(reslow)) { 62 stop("reslow is not a vector") 63 } 64 if(length(reslow) != k) { 65 stop("reslow has not the right dimension") 66 } 67 } 68 if(!is.null(reshigh)) { 69 if(!is.vector(reshigh)) { 70 stop("reshigh is not a vector") 71 } 72 if(length(reshigh) != k) { 73 stop("reshigh has not the right dimension") 74 } 75 } 76 if(riskless) { 77 a1 <- colMeans(x) - rf 78 if(shorts) { 79 a2 <- NULL 80 b2 <- NULL 81 } 82 else { 83 a2 <- matrix(0, k, k) 84 diag(a2) <- 1 85 b2 <- rep.int(0, k) 86 } 87 if(!is.null(reslow) && !is.null(reshigh)) { 88 a3 <- matrix(0, k, k) 89 diag(a3) <- 1 90 Amat <- t(rbind(a1, a2, a3, -a3)) 91 b0 <- c(pm-rf, b2, reslow, -reshigh) 92 } 93 else { 94 Amat <- t(rbind(a1, a2)) 95 b0 <- c(pm-rf, b2) 96 } 97 res <- solve.QP(Dmat, dvec, Amat, bvec=b0, meq=1) 98 } 99 else { 100 a1 <- rep.int(1, k) 101 a2 <- colMeans(x) 102 if(shorts) { 103 if(!is.null(reslow) && !is.null(reshigh)) { 104 a3 <- matrix(0, k, k) 105 diag(a3) <- 1 106 Amat <- t(rbind(a1, a2, a3, -a3)) 107 b0 <- c(1, pm, reslow, -reshigh) 108 } 109 else { 110 Amat <- t(rbind(a1, a2)) 111 b0 <- c(1, pm) 112 } 113 } 114 else { 115 a3 <- matrix(0, k, k) 116 diag(a3) <- 1 117 b3 <- rep.int(0, k) 118 if(!is.null(reslow) && !is.null(reshigh)) { 119 Amat <- t(rbind(a1, a2, a3, a3, -a3)) 120 b0 <- c(1, pm, b3, reslow, -reshigh) 121 } 122 else { 123 Amat <- t(rbind(a1, a2, a3)) 124 b0 <- c(1, pm, b3) 125 } 126 } 127 res <- solve.QP(Dmat, dvec, Amat, bvec=b0, meq=2) 128 } 129 y <- c(tcrossprod(res$solution, x)) 130 ans <- list(pw = res$solution, px = y, pm = mean(y), ps = sd(y)) 131 return(ans) 132} 133 134get.hist.quote <- 135function(instrument = "^gdax", start, end, 136 quote = c("Open", "High", "Low", "Close"), 137 provider = c("yahoo"), method = NULL, 138 origin = "1899-12-30", compression = "d", 139 retclass = c("zoo", "ts"), 140 quiet = FALSE, drop = FALSE) 141{ 142 if(missing(start)) start <- "1991-01-02" 143 if(missing(end)) end <- format(Sys.Date() - 1, "%Y-%m-%d") 144 145 provider <- match.arg(provider) 146 retclass <- match.arg(retclass) 147 148 periodicity <- match.arg(compression, 149 c("daily", "weekly", "monthly")) 150 151 ## Be nice. 152 ind <- pmatch(quote, "AdjClose", nomatch = 0L) 153 quote[ind] <- "Adjusted" 154 155 start <- as.Date(start) 156 end <- as.Date(end) 157 158 x <- getSymbols(instrument, src = "yahoo", 159 from = start, to = end, return.class = "zoo", 160 periodicity = periodicity, auto.assign = FALSE) 161 162 colnames(x) <- sub(".*\\.", "", colnames(x)) 163 nser <- pmatch(quote, colnames(x)) 164 if(any(is.na(nser))) 165 stop("this quote is not available") 166 n <- nrow(x) 167 168 dat <- index(x) 169 if(!quiet && (dat[1] != start)) 170 cat(format(dat[1], "time series starts %Y-%m-%d\n")) 171 if(!quiet && (dat[n] != end)) 172 cat(format(dat[n], "time series ends %Y-%m-%d\n")) 173 174 if(retclass == "ts") { 175 jdat <- unclass(julian(dat, origin = as.Date(origin))) 176 ## We need unclass() because 1.7.0 does not allow adding a 177 ## number to a "difftime" object. 178 ind <- jdat - jdat[1] + 1 179 y <- matrix(NA, nrow = max(ind), ncol = length(nser)) 180 y[ind, ] <- as.matrix(x[, nser, drop = FALSE]) 181 colnames(y) <- colnames(x)[nser] 182 y <- y[, seq_along(nser), drop = drop] 183 return(ts(y, start = jdat[1], end = jdat[n])) 184 } else { 185 x[ , nser, drop = drop] 186 } 187} 188 189maxdrawdown <- 190function(x) 191{ 192 if(NCOL(x) > 1) 193 stop("x is not a vector or univariate time series") 194 if(any(is.na(x))) 195 stop("NAs in x") 196 cmaxx <- cummax(x)-x 197 mdd <- max(cmaxx) 198 to <- which(mdd == cmaxx) 199 from <- double(NROW(to)) 200 for (i in 1:NROW(to)) 201 from[i] <- max(which(cmaxx[1:to[i]] == 0)) 202 return(list(maxdrawdown = mdd, from = from, to = to)) 203} 204 205sharpe <- 206function(x, r = 0, scale = sqrt(250)) 207{ 208 if(NCOL(x) > 1) 209 stop("x is not a vector or univariate time series") 210 if(any(is.na(x))) 211 stop("NAs in x") 212 if(NROW(x) == 1) 213 return(NA) 214 else { 215 y <- diff(x) 216 return(scale * (mean(y)-r)/sd(y)) 217 } 218} 219 220sterling <- 221function(x) 222{ 223 if(NCOL(x) > 1) 224 stop("x is not a vector or univariate time series") 225 if(any(is.na(x))) 226 stop("NAs in x") 227 if(NROW(x) == 1) 228 return(NA) 229 else { 230 return((x[NROW(x)]-x[1]) / maxdrawdown(x)$maxdrawdown) 231 } 232} 233 234plotOHLC <- 235function(x, xlim = NULL, ylim = NULL, xlab = "Time", ylab, 236 col = par("col"), bg = par("bg"), axes = TRUE, 237 frame.plot = axes, ann = par("ann"), main = NULL, 238 date = c("calendar", "julian"), format = "%Y-%m-%d", 239 origin = "1899-12-30", ...) 240{ 241 if ((!is.mts(x)) || 242 (colnames(x)[1] != "Open") || 243 (colnames(x)[2] != "High") || 244 (colnames(x)[3] != "Low") || 245 (colnames(x)[4] != "Close")) 246 stop("x is not a open/high/low/close time series") 247 xlabel <- if (!missing(x)) 248 deparse(substitute(x)) 249 else NULL 250 if (missing(ylab)) 251 ylab <- xlabel 252 date <- match.arg(date) 253 time.x <- time(x) 254 dt <- min(lag(time.x)-time.x)/3 255 if (is.null(xlim)) 256 xlim <- range(time.x) 257 if (is.null(ylim)) 258 ylim <- range(x[is.finite(x)]) 259 plot.new() 260 plot.window(xlim, ylim, ...) 261 segments(time.x, x[, "High"], time.x, x[, "Low"], col = col[1], bg = bg) 262 segments(time.x - dt, x[, "Open"], time.x, x[, "Open"], 263 col = col[1], bg = bg) 264 segments(time.x, x[, "Close"], time.x + dt, x[, "Close"], 265 col = col[1], bg = bg) 266 if (ann) 267 title(main = main, xlab = xlab, ylab = ylab, ...) 268 if (axes) { 269 if (date == "julian") { 270 axis(1, ...) 271 axis(2, ...) 272 } 273 else { 274 n <- NROW(x) 275 lab.ind <- round(seq(1, n, length.out = 5)) 276 D <- as.vector(time.x[lab.ind]*86400) + as.POSIXct(origin, tz = "GMT") 277 DD <- format.POSIXct(D, format = format, tz ="GMT") 278 axis(1, at=time.x[lab.ind], labels=DD, ...) 279 axis(2, ...) 280 } 281 } 282 if (frame.plot) 283 box(...) 284} 285 286