pp.convdate <- function (x) { x <- as.character(x) xs <- strsplit(x, " ") days <- sapply(xs, "[", 1) mons <- sapply(xs, "[", 2) yrs <- sapply(xs, "[", 3) mons <- match(mons, c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December")) mons <- as.character(mons) nm <- nchar(mons) mons[nm == 1] <- paste("0", mons[nm == 1], sep = "") nd <- nchar(days) days[nd == 1] <- paste("0", days[nd == 1], sep = "") paste(yrs, mons, days, sep = "-") } pp.denplot <- function (ret, name, start, end=NULL, lwd=2, main="", tlim=3, probs=c(.99, .95, .90, .75, .25, .10, .05, .01), digits=0, xlim=range(x, start, end), ...) { # function placed in the public domain 2011 by Burns Statistics x <- start * exp(ret) denx <- density(x, to=min(max(x), tlim * start)) plot(denx, main=main, xlab=name, ylab="", lwd=lwd, col="steelblue", yaxt="n", xlim=xlim, ...) abline(v=start, col="green") if(length(end)) abline(v=end, col="red") quans <- quantile(x, probs=probs) qtext <- paste(format(probs*100), "%: ", format(round(quans, digits=digits)), sep="") text(max(denx$x) - .02 * start, adj=1, .8 * par("usr")[4], paste(qtext, collapse="\n")) } pp.garchblocksim <- function (garobj, block, ahead = 250, number = 5000, collapse = TRUE, var.out = FALSE, trim = 0.001) { subfun.samp <- function(resids, top, bseq, nb, ahead, block) { ans <- rep(NA * 0, ahead) for (i in 1:nb) { this.samp <- resids[sample(top, 1) + bseq] if (runif(1) < 0.5) this.samp <- -this.samp ans[bseq + 1 + block * (i - 1)] <- this.samp } length(ans) <- ahead ans } if (length(garobj$coef) != 3) stop("garch(1,1) model only") resids <- garobj$residuals resids <- resids[!is.na(resids)] top <- length(resids) - block + 1 bseq <- 0:(block - 1) nb <- ceiling(ahead/block) coefs <- garobj$coef int <- coefs[1] alph <- coefs[2] beta <- coefs[3] startvar <- garobj$fitted.values[length(resids)]^2 out <- array(NA, c(ahead, number)) if (var.out) outv <- out for (j in 1:number) { this.res <- subfun.samp(resids, top, bseq, nb, ahead, block) this.var <- startvar for (i in 1:ahead) { out[i, j] <- this.ret <- sqrt(this.var) * this.res[i] this.var <- int + alph * this.ret^2 + beta * this.var if (var.out) outv[i, j] <- this.var } } if (collapse) { out <- colSums(out) if (trim > 0) { tlim <- quantile(out, c(trim, 1 - trim)) out <- out[out > tlim[1] & out < tlim[2]] } } if (var.out) { list(returns = out, condvar = outv, garch.coef = coefs) } else { attr(out, "garch.coef") <- coefs out } } pp.garchsim <- function (garobj, ahead = 250, number = 5000, collapse = TRUE, var.out = FALSE, trim = 0.001) { if (length(garobj$coef) != 3) stop("garch(1,1) model only") resids <- garobj$residuals resids <- resids[!is.na(resids)] resids <- c(resids, -resids) coefs <- garobj$coef int <- coefs[1] alph <- coefs[2] beta <- coefs[3] startvar <- garobj$fitted.values[length(resids)]^2 out <- array(NA, c(ahead, number)) if (var.out) outv <- out for (j in 1:number) { this.res <- sample(resids, ahead, replace = TRUE) this.var <- startvar for (i in 1:ahead) { out[i, j] <- this.ret <- sqrt(this.var) * this.res[i] this.var <- int + alph * this.ret^2 + beta * this.var if (var.out) outv[i, j] <- this.var } } if (collapse) { out <- colSums(out) if (trim > 0) { tlim <- quantile(out, c(trim, 1 - trim)) out <- out[out > tlim[1] & out < tlim[2]] } } if (var.out) { list(returns = out, condvar = outv, garch.coef = coefs) } else { attr(out, "garch.coef") <- coefs out } } pp.htmltable <- function (x, filename = "", order = TRUE) { sfun.row <- function(r) { cat(file = filename, append = TRUE, "\n") cat(file = filename, append = TRUE, paste("", r, "\n", sep = "")) cat(file = filename, append = TRUE, "\n") } if (!is.data.frame(x)) stop("expecting 'x' to be a data frame") if (ncol(x) != 4) stop("expecting 'x' to have 4 columns") if (order) { x <- x[order(x[, "value"]), ] } ref <- as.character(x[, "reference"]) ref <- paste("webpage", sep = "") x[, "reference"] <- ref x[, "value"] <- format(x[, "value"], big.mark = ",") x <- as.matrix(x) cat("\n\n", file = filename) sfun.row(paste("", colnames(x), "", sep = "")) for (i in 1:nrow(x)) sfun.row(x[i, ]) cat("\n
\n", file = filename, append = TRUE) } pp.lores <- function (x, nobs) { n <- length(x) if (nobs > n) { res <- x - mean(x) } else { df <- data.frame(time = 1:n, y = drop(x)) lo <- loess(y ~ time, data = df, span = nobs/n) res <- lo$residuals } res } pp.omnipred <- function (symbol, name, prefix, year=2012, startdate=20030101, nobs=7, block=50, lev=NULL, garsave=FALSE, trim=TRUE, Debug=FALSE) { # function placed in the public domain 2011 by Burns Statistics require(TTR) require(tseries) if(!length(lev)) { lev <- getYahooData(symbol, start=startdate)[, "Close"] } if(trim) { dat <- index(lev) out <- substring(dat, 1, 4) == year lev <- lev[!out] } ret <- diff(log(lev)) res <- pp.lores(ret, nobs) sco1 <- c(var(res) * .02, .05, .93) gar1 <- garch(res, start=sco1) sco2 <- c(var(res) * .01, .05, .94) gar2 <- garch(res, start=sco2) gsum1 <- sum(gar1$coef[-1]) gsum2 <- sum(gar2$coef[-1]) cat("fraction of zeros in res is", mean(abs(unclass(res)) < 1e-10), "\n") cat("garch coeficients 1 are:\n") print(gar1$coef) cat("alpha + beta 1:", gsum1, "\n") cat("garch coeficients 2 are:\n") print(gar2$coef) cat("alpha + beta 2:", gsum2, "\n") initlev <- as.vector(tail(lev, 1)) cat("initial level is", initlev, "\n") if(gsum1 < 1 && gsum2 >= 1) { gar <- gar1 } else if(gsum2 < 1 && gsum1 >= 1) { gar <- gar2 } else { if(gar1$n.likeli < gar2$n.likeli) { gar <- gar1 } else { gar <- gar2 } } fullname <- paste(name, year, sep=" -- ") if(block > 1) { sim <- pp.garchblocksim(gar, block) } else { sim <- pp.garchsim(gar) } if(garsave) { assign(paste(prefix, year, "gar", sep="_"), gar, pos=1) } assign(paste(prefix, year, "sim", sep="_"), sim, pos=1) assign(paste(prefix, year, "res", sep="_"), res, pos=1) if(Debug) assign(paste(prefix, year, "lev", sep="_"), lev, pos=1) assign(paste(prefix, year, "initlev", sep="_"), initlev, pos=1) assign(paste(prefix, year, "ecdf", sep="_"), ecdf(initlev * exp(sim)), pos=1) plotfile <- paste(prefix, "_", year, ".png", sep="") png(file=plotfile, width=512) on.exit(dev.off()) par(mar=c(4,2,0,2) + .1) pp.denplot(sim, fullname, start=initlev) plotfile } pp.png <- function (filename, mar = c(5, 4, 0, 2), ...) { png(file = filename, width = 512) par(mar = mar + 0.1, ...) } pp.predframe <- function (prefix, person, firm, value, reference, note = "", year = 2011) { pfnam <- paste(prefix, "_", year, ".predframe", sep = "") nin <- length(person) if (length(firm) != nin) stop("bad length for 'firm'") if (length(value) != nin) stop("bad length for 'value'") if (length(reference) != nin) stop("bad length for 'reference'") note <- rep(note, length = nin) new <- data.frame(person = person, firm = firm, value = value, reference = reference) if (exists(pfnam, where = 1)) { old <- get(pfnam, pos = 1) ans <- rbind(old, new) attr(ans, "prefix") <- attr(old, "prefix") attr(ans, "year") <- attr(old, "year") attr(ans, "description") <- attr(old, "description") } else { ans <- new attr(ans, "prefix") <- prefix attr(ans, "year") <- year cat("Need to add 'description' attribute to", pfnam, "\n") } assign(pfnam, ans, pos = 1) } pp.predplot <- function (pframe, ..., end=NULL) { # function placed in the public domain 2011-2012 by Burns Statistics fprefix <- paste(attr(pframe, "prefix"), attr(pframe, "year"), sep="_") name <- paste(attr(pframe, "description"), attr(pframe, "year"), sep=" - ") ret <- get(paste(fprefix, "sim", sep="_")) start <- get(paste(fprefix, "initlev", sep="_")) lrange <- start * exp(range(ret)) pp.denplot(ret, name, start, end=end, xlim=range(lrange, pframe[, "value"], end), ...) abline(v=pframe[, "value"], col="purple") } pp.yearendplot <- function (name, prefix, year=2011) { # function placed in the public domain 2012 by Burns Statistics fullname <- paste(name, year, sep=" -- ") sim <- get(paste(prefix, year, "sim", sep="_")) initlev <- get(paste(prefix, year, "initlev", sep="_")) endlev <- get(paste(prefix, year + 1, "initlev", sep="_")) pfn <- paste(prefix, '_', year, '.predframe', sep='') plotfile <- paste(prefix, "_", year, "_end.png", sep="") png(file=plotfile, width=512) on.exit(dev.off()) par(mar=c(4,2,0,2) + .1) if(exists(pfn, where=1)) { pp.predplot(get(pfn), end=endlev) } else { pp.denplot(sim, fullname, start=initlev, end=endlev) } plotfile }