[Uwgarp-commits] r99 - pkg/GARPFRM/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Mar 5 04:23:12 CET 2014
Author: rossbennett34
Date: 2014-03-05 04:23:09 +0100 (Wed, 05 Mar 2014)
New Revision: 99
Modified:
pkg/GARPFRM/R/EWMA.R
pkg/GARPFRM/R/monte_carlo.R
Log:
Modifying EWMA source code and documentation
Modified: pkg/GARPFRM/R/EWMA.R
===================================================================
--- pkg/GARPFRM/R/EWMA.R 2014-03-05 02:23:10 UTC (rev 98)
+++ pkg/GARPFRM/R/EWMA.R 2014-03-05 03:23:09 UTC (rev 99)
@@ -1,51 +1,19 @@
-# EWMA <- function(object, lambda=0.96, cor=FALSE) {
-# if ((lambda<1 || lambda > 0)){
-# object.names = colnames(object)
-# t.object = nrow(object)
-# k.object = ncol(object)
-# object = as.matrix(object)
-# t.names = rownames(object)
-#
-# covEWMA = array(,c(t.object,k.object,k.object))
-# # Note it is unconditional cov
-# cov.f = var(object)
-# FF = (object[1,]- mean(object)) %*% t(object[1,]- mean(object))
-# covEWMA[1,,] = (1-lambda)*FF + lambda*cov.f
-# for (i in 2:t.object) {
-# FF = (object[i,]- mean(object)) %*% t(object[i,]- mean(object))
-# print(FF)
-# covEWMA[i,,] = (1-lambda)*FF + lambda*covEWMA[(i-1),,]
-# }
-#
-# } else {
-# stop("exp-decay lambda must be ]0:1[")
-# }
-#
-# dimnames(covEWMA) = list(t.names, object.names, object.names)
-#
-# if(cor) {
-# corEWMA = covEWMA
-# for (i in 1:dim(corEWMA)[1]) {
-# corEWMA[i, , ] = cov2cor(covEWMA[i, ,])
-# }
-# return(corEWMA)
-# } else{
-# return(covEWMA)
-# }
-# }
+
#' Exponential Weighted Moving Average (EWMA)
#'
-#' Description of EWMA. The function handles UV and MLM objects and returns either cov/cor.
+#' Use an exponentially weighted moving average to estimate the covariance or
+#' correlation of asset returns.
#'
-#' @param R
-#' @param lambda
-#' @param initialWindow is the initializing window
-#' @param correlation option (cor by default = FALSE)
+#' @param R asset returns
+#' @param lambda smoothing parameter, must be greater than 0 or less than 1
+#' @param initialWindow initial window of observations used in estimating the
+#' initial covariance or correlation
+#' @param TRUE/FALSE to return a correlation matrix. Default cor = FALSE.
#' @export
EWMA <- function(R, lambda=0.94, initialWindow=10, cor=FALSE){
# Check for lambda between 0 and 1 & initialWindow must be greater than ncol(R)
- if (((lambda>1 || lambda < 0))) {stop("For exponential decay lambda must belong to ]0:1[")}
- if (initialWindow> nrow(R)){stop("Initialization window is too large")}
+ if ((lambda > 1) || (lambda < 0)) stop("lambda must be greater than 0 or less than 1")
+ if (initialWindow > nrow(R)) stop("Initial window must be less than the number of observations")
# Separate data into a initializing window and a testing window
initialR = R[1:initialWindow,]
@@ -73,111 +41,101 @@
if(cor & ncol(R) > 1) {
out$EWMA <- lapply(out$EWMA, cov2cor)
class(out) <- c("EWMACor")
- } else if(cor & ncol(R)==1) {
+ } else if(cor & (ncol(R) == 1)) {
stop("EWMA correlation is only to be estimated for two or more assets")
}
# Check for Covar or Var
if((cor == FALSE) & (ncol(R) > 1)) {
- class(out) <- c("EWMACovar")
+ class(out) <- c("covEWMA")
} else if ((cor == FALSE) & (ncol(R) == 1)){
- class(out) <- c("EWMAVar")
+ class(out) <- c("varEWMA")
}
- # Bind initial data to EWMA object in order to plot a comparison
- # out$y_data <- R
- # out$y_data <- R adds R to the list element
- # The EWMA estimate and R should be separate elements in the list returned
return(out)
}
-#' EWMA Volatility/Cross-Volatility
+#' EWMA Covariance
#'
-#' Description of EWMA Vola
+#' Extract the covariance between two assets from an EWMA object
#'
-#' @param object a EWMA object created by \code{\link{EWMA}}
+#' @param EWMA an EWMA object created by \code{\link{EWMA}}
+#' @param assets character vector or numeric vector. If
+#' \code{assets} is of length 1, then the variance will be returned.
+#' The assets can be specified by name or index.
#' @export
-getCov <- function(object, asset1, asset2){
+getCov <- function(EWMA, assets){
UseMethod("getCov")
}
-#' @method getCov EWMACovar
-#' @S3method getCov EWMACovar
-getCov.EWMACovar <- function(object, asset1, asset2){
- if(!inherits(object, "EWMACovar")) stop("object must be of class EWMACovar")
- # Manipulate object for feasible use
- # object[[length(object)]] = NULL
+#' @method getCov covEWMA
+#' @S3method getCov covEWMA
+getCov.covEWMA <- function(EWMA, assets=c(1,2)){
+ if(!inherits(EWMA, "covEWMA")) stop("object must be of class covEWMA")
# Check if asset is a character
- if(is.character(asset1) & is.character(asset2)){
- idx1 = grep(asset1, colnames(object$EWMA[[1]]))
+ if(is.character(assets[1]) & is.character(assets[2])){
+ idx1 = grep(assets[1], colnames(EWMA$EWMA[[1]]))
if(length(idx1) == 0) stop("name for asset1 not in object")
- idx2 = grep(asset2, colnames(object$EWMA[[1]]))
+ idx2 = grep(assets[2], colnames(EWMA$EWMA[[1]]))
if(length(idx2) == 0) stop("name for asset2 not in object")
} else {
# Then dimensions are enough to find covar
- idx1 = asset1
- idx2 = asset2
+ idx1 = assets[1]
+ idx2 = assets[1]
}
- out = xts(unlist(lapply(object$EWMA, function(X) X[idx1, idx2])), as.Date(names(object$EWMA)))
- colnames(out) = paste(asset1, asset2, sep=".")
+ out = xts(unlist(lapply(EWMA$EWMA, function(X) X[idx1, idx2])), as.Date(names(EWMA$EWMA)))
+ colnames(out) = paste(assets[1], assets[2], sep=".")
return(out)
}
-#' @method getCov EWMAVar
-#' @S3method getCov EWMAVar
-getCov.EWMAVar <- function(object, asset1, asset2=NULL){
- if(!inherits(object, "EWMAVar")) stop("Object must be of class EWMAVar")
- if (is.null(asset2) == FALSE) {warning("Running univariate no need to specify asset2")}
- # Manipulate object for feasible use
- # object[[length(object)]] = NULL
+#' @method getCov varEWMA
+#' @S3method getCov varEWMA
+getCov.varEWMA <- function(EWMA, assets=1){
+ if(!inherits(EWMA, "varEWMA")) stop("EWMA must be of class varEWMA")
# Check if asset is a character
- if(is.character(asset1)){
- idx1 = grep(asset1, colnames(object$EWMA[[1]]))
- if(length(idx1) == 0) stop("name for asset1 not in object")
+ if(is.character(assets[1])){
+ idx1 = grep(assets[1], colnames(EWMA$EWMA[[1]]))
+ if(length(idx1) == 0) stop("name for asset not in EWMA object")
} else {
- # Then dimensions are enough to find covar
- idx1 = asset1
+ idx1 = assets[1]
}
- out = xts(unlist(lapply(object$EWMA, function(x) x[idx1])), as.Date(names(object$EWMA)))
- colnames(out) = asset1
-
+ out = xts(unlist(lapply(EWMA$EWMA, function(x) x[idx1])), as.Date(names(EWMA$EWMA)))
+ colnames(out) = assets[1]
return(out)
}
#' EWMA Correlation
#'
-#' Description of EWMA Correlation, requires two assets
+#' Extract the correlation of two assets from an EWMA object
#'
-#' @param object a EWMA object created by \code{\link{EWMA}}
+#' @param object an EWMA object created by \code{\link{EWMA}}
#' @export
-getCor <- function(object, asset1, asset2){
+getCor <- function(EWMA, assets){
UseMethod("getCor")
}
-#' @method getCor EWMACor
-#' @S3method getCor EWMACor
-getCor.EWMACor <- function(object, asset1, asset2){
- if(!inherits(object, "EWMACor")) stop("object must be of class EWMACor")
- # Manipulate object for feasible use
- # object[[length(object)]] = NULL
+#' @method getCor corEWMA
+#' @S3method getCor corEWMA
+getCor.corEWMA <- function(EWMA, assets=c(1,2)){
+ if(!inherits(EWMA, "corEWMA")) stop("EWMA must be of class corEWMA")
# Check if asset is a character
- if(is.character(asset1) & is.character(asset2)){
- idx1 = grep(asset1, colnames(object$EWMA[[1]]))
- if(length(idx1) == 0) stop("name for asset1 not in object")
- idx2 = grep(asset2, colnames(object$EWMA[[1]]))
- if(length(idx2) == 0) stop("name for asset2 not in object")
+ if(is.character(assets[1]) & is.character(assets[2])){
+ idx1 = grep(assets[1], colnames(EWMA$EWMA[[1]]))
+ if(length(idx1) == 0) stop("name for asset1 not in EWMA object")
+ idx2 = grep(assets[2], colnames(EWMA$EWMA[[1]]))
+ if(length(idx2) == 0) stop("name for asset2 not in EWMA object")
} else {
# Then dimensions are enough to find covar
- idx1 = asset1
- idx2 = asset2
+ idx1 = assets[1]
+ idx2 = assets[2]
}
- out = xts(unlist(lapply(object$EWMA, function(x) x[idx1, idx2])), as.Date(names(object$EWMA)))
- colnames(out) = paste(asset1, asset2, sep=".")
+ out = xts(unlist(lapply(EWMA$EWMA, function(x) x[idx1, idx2])), as.Date(names(EWMA$EWMA)))
+ colnames(out) = paste(assets[1], assets[2], sep=".")
return(out)
}
@@ -190,21 +148,21 @@
# EWMA plotting for covar
#' @export
-plot.EWMACovar <- function(object,..., asset1, asset2){
+plot.covEWMA <- function(object, ..., assets=c(1, 2)){
# Check if asset is a character
- if(is.character(asset1) & is.character(asset2)){
- idx1 = grep(asset1, colnames(object[[1]]))
- if(length(idx1) == 0) stop("name for asset1 not in object")
- idx2 = grep(asset2, colnames(object[[1]]))
- if(length(idx2) == 0) stop("name for asset2 not in object")
+ if(is.character(assets[1]) & is.character(assets[2])){
+ idx1 = grep(assets[1], colnames(object[[1]]))
+ if(length(idx1) == 0) stop("name for assets[1] not in object")
+ idx2 = grep(assets[2], colnames(object[[1]]))
+ if(length(idx2) == 0) stop("name for assets[2] not in object")
} else {
# Then dimensions are enough to find covar
- idx1 = asset1
- idx2 = asset2
+ idx1 = assets[1]
+ idx2 = assets[2]
}
- tmp = getCov(object,..., asset1, asset2)
- plot(x=time(as.zoo(tmp)), y=tmp, type="l", xlab="Time", ylab="Covariance", lwd=2, col="blue",
- main="EWMA Covariance");
+ tmp = getCov(object, ..., assets)
+ plot(x=time(as.zoo(tmp)), y=tmp, type="l", xlab="Time", ylab="Covariance",
+ lwd=2, col="blue", main="EWMA Covariance")
grid()
abline(h=var(object$R)[idx1,idx2], lwd=2, col="red")
}
@@ -212,10 +170,10 @@
# EWMA plotting for var
#' @export
-plot.EWMAVar <- function(object,...,asset1){
- tmp = getCov(object,asset1)
- plot(x=time(as.zoo(tmp)),y=tmp, type="l", xlab="Time", ylab="Variance", lwd=2, col="blue",
- main="EWMA Variance");
+plot.varEWMA <- function(object, ..., assets=c(1,2)){
+ tmp = getCov(object, assets[1])
+ plot(x=time(as.zoo(tmp)),y=tmp, type="l", xlab="Time", ylab="Variance",
+ lwd=2, col="blue", main="EWMA Variance");
grid()
abline(h=var(object$R), lwd=2, col="red")
}
@@ -223,21 +181,21 @@
# EWMA plotting for correlation
#' @export
-plot.EWMACor <- function(object,...,asset1, asset2){
+plot.corEWMA <- function(object, ..., assets=c(1,2)){
# Check if asset is a character
- if(is.character(asset1) & is.character(asset2)){
- idx1 = grep(asset1, colnames(object[[1]]))
+ if(is.character(assets[1]) & is.character(assets[2])){
+ idx1 = grep(assets[1], colnames(object[[1]]))
if(length(idx1) == 0) stop("name for asset1 not in object")
- idx2 = grep(asset2, colnames(object[[1]]))
+ idx2 = grep(assets[2], colnames(object[[1]]))
if(length(idx2) == 0) stop("name for asset2 not in object")
} else {
# Then dimensions are enough to find covar
- idx1 = asset1
- idx2 = asset2
+ idx1 = assets[1]
+ idx2 = assets[2]
}
- tmp = getCor(object, asset1, asset2)
- plot(x=time(as.zoo(tmp)), y=tmp, type="l", xlab="Time", ylab="Correlation", lwd=2, col="blue",
- main="EWMA Correlation");
+ tmp = getCor(object, assets)
+ plot(x=time(as.zoo(tmp)), y=tmp, type="l", xlab="Time", ylab="Correlation",
+ lwd=2, col="blue", main="EWMA Correlation")
grid()
abline(h=cor(object$R)[idx1,idx2], lwd=2, col="red")
}
Modified: pkg/GARPFRM/R/monte_carlo.R
===================================================================
--- pkg/GARPFRM/R/monte_carlo.R 2014-03-05 02:23:10 UTC (rev 98)
+++ pkg/GARPFRM/R/monte_carlo.R 2014-03-05 03:23:09 UTC (rev 99)
@@ -1,10 +1,11 @@
# Monte Carlo Function
generateMC <- function(mu, sigma, Time=1, steps=52, starting_value=100){
dt <- Time / steps
- S <- vector("numeric", steps)
+ S <- vector("numeric", steps+1)
+ eps <- rnorm(steps)
S[1] <- starting_value
for(i in 2:length(S)){
- dS <- mu * dt + sigma * rnorm(1) * sqrt(dt)
+ dS <- mu * dt + sigma * eps(i-1) * sqrt(dt)
S[i] <- dS + S[i-1]
}
return(S)
@@ -14,27 +15,38 @@
# more accurate
generateLogMC <- function(mu, sigma, Time=1, steps=52, starting_value=100){
dt <- Time / steps
- S <- vector("numeric", steps)
+ S <- vector("numeric", steps+1)
S[1] <- starting_value
for(i in 2:length(S)){
- S[i] <- S[i-1] * exp((mu - 0.5 * sigma^2) * dt + sigma * rnorm(1) * sqrt(dt))
+ S[i] <- S[i-1] * exp((mu - 0.5 * sigma^2) * dt + sigma * eps(i-1) * sqrt(dt))
}
return(S)
}
#' Monte Carlo Price Path Simulation
#'
-#' Description for Monte Carlo. Geometric brownian motion. mu and sigma are assumed constant
+#' Run \code{N} monte carlo simulations to generate asset price paths following
+#' a geometric brownian motion process.
#'
+#' TODO: add equations for GBM
+#'
+#' @note This function returns a m x N matrix of simulated price paths where
+#' m is the number of steps + 1 and N is the number of simulations. This can be
+#' very memory and compute intensive with a large number of steps and/or a
+#' large number of simulations.
+#' More efficient methods in terms of speed and memory should be used, for
+#' example, to price options.
+#'
#' @param mu annualized expected return
#' @param sigma annualized standard deviation
#' @param N number of simulations
#' @param Time length of simulation (in years)
#' @param steps number of time steps
-#' @param starting_value price to start at
-#' @param log TRUE/FALSE (default = TRUE) simulate ln(P) rather than S; where S
+#' @param starting_value asset price starting value
+#' @param log TRUE/FALSE (default = TRUE) simulate ln(S) rather than S; where S
#' is the price of the asset.
-#' @return matrix of Monte Carlo simulated price paths
+#' @return matrix of simulated price paths where each column represents a price path
+#' @export
monteCarlo <- function(mu, sigma, N=100, Time=1, steps=52, starting_value=100, log=TRUE){
mc_mat <- matrix(0, steps, N)
if(log){
@@ -51,7 +63,7 @@
}
plot.monte_carlo <-function(x, y, ..., main="Monte Carlo Simulation", xlab="Time Index", ylab="Price"){
- plot(x[,1], type="n", ylim=range(x), main=main, xlab=xlab, ylab=ylab)
+ plot(x[,1], type="n", ylim=range(x), main=main, xlab=xlab, ylab=ylab, ...)
for(i in 1:ncol(x)){
lines(x[,i])
}
@@ -59,11 +71,22 @@
#' Ending Prices of Monte Carlo Simulation
#'
+#' Get the ending prices, i.e. terminal values, of a monte carlo simulation
+#' @param mc monte carlo object created with \link{\code{monteCarlo}}
+#' @return vector ending prices
+#' @export
endingPrices <- function(mc){
- mc[nrow(mc),]
+ if(!inherits(mc, "monte_carlo")) stop("mc must be of class 'monte_carlo'")
+ return(mc[nrow(mc),])
}
+#' Plot Ending Prices
+#'
+#' Plot the ending prices from a Monte Carlo simulation
+#' @param mc monte carlo object created with \link{\code{monteCarlo}}
+#' @export
plotEndingPrices <- function(mc){
+ if(!inherits(mc, "monte_carlo")) stop("mc must be of class 'monte_carlo'")
ending_prices <- endingPrices(mc)
dens_ep <- density(ending_prices)
hist(ending_prices, freq=FALSE)
More information about the Uwgarp-commits
mailing list