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