# Load the required packages. | |
library(PerformanceAnalytics) | |
library(xts) | |
#' @title Sample a random disturbance from either a normal distribution with a | |
#' constant standard deviation (Geometric Brownian Motion model) or from a | |
#' distribution with a stochastic standard deviation (Stochastic Volatility GBM) | |
#' | |
#' @description Given a long run random disturbance mean, mu, and a standard | |
#' deviation, sigma, sample a random disturbance from either a standard normal | |
#' distribution with mean mu and standard deviation sigma, or sample a random | |
#' disturbance from a normal distribution with mean mu and a stochastic standard | |
#' deviation sampled from a normal distribution centered around sigma with a | |
#' standard deviation equal to half of sigma. This is maxed out at 0.01 because | |
#' a negative standard deviation is incoherent. | |
#' | |
#' @details When the stochastic volatility option is enabled, the resulting | |
#' distribution of means for the random disturbance has much fatter tails. This | |
#' can be illustrated with the compareRandomDisturbances function. | |
#' | |
#' @param mu double :: The expected value of the disturbance (should be zero). | |
#' @param sigma double :: The standard deviation of the disturbance. | |
#' | |
randomDisturbance <- function(mu = 0.0, sigma = 0.45, | |
stochastic.volatility = TRUE) { | |
if (!stochastic.volatility) { | |
return (rnorm(1, mean = mu, sd = sigma)) | |
} else { | |
# sigma.stochastic <- rlnorm(1, meanlog = log(sigma), sdlog = sigma/2) | |
sigma.stochastic <- max(rnorm(1, mean = sigma, sd = sigma/2), 0.01) | |
return (rnorm(1, mean = mu, sd = sigma.stochastic)) | |
} | |
} | |
#' @title Step forward one unit of time in the log price process. | |
#' @description This step function is used to move from state t / t-1 in the log | |
#' price process to t+1 / t. This involves adding a static drift component to the | |
#' previous price and adding a random disturbance sampled from some distribution. | |
#' | |
#' @param Xt1 double :: The current price in the log price process. | |
#' @param mu double :: The drift component of returns. | |
#' @param rd double :: A random disturbance in the force. | |
#' @param dt double :: This variable controls time. It is set to 1/252 by default | |
#' @return Xt2 double :: The next price in the log price process. | |
#' | |
logPriceStep <- function(Xt1, mu, rd, dt = 0.003968254) { | |
Xt2 <- Xt1 + (mu * dt) + (rd * sqrt(dt)) | |
return (Xt2) | |
} | |
#' @title Generate a log price process. | |
#' @description This method is used to generate a log price process which | |
#' simulates the log price of an asset starting with the value X0. The process | |
#' is controlled by a number of parameters including mu, the average daily | |
#' return, rd.mu the average daily random disturbance, rd.sigma, the volatility | |
#' of the random disturbance, and dt, time. The method works by iteratively | |
#' sampling a new log price using the logPriceStep function. | |
#' | |
#' @param t int :: The number of days worth of prices to simulate. | |
#' @param X0 int :: The starting price of the asset. Set to 1.0. | |
#' @param mu double :: The drift component of returns. | |
#' @param rd.mu double :: The average daily disturbance to be added. | |
#' @param rd.sigma double :: The volatility of daily disturbances to be added. | |
#' @param dt double :: Time. This variable is set to 1/252 by default. | |
#' @return X vector :: A simulated log price process. | |
#' | |
logPriceProcess <- function(t = 252, X0 = 1.0, mu = 0.1062, | |
rd.mu = 0.0, rd.sigma = 0.45, | |
dt = 0.003968254, stochastic.volatility = TRUE) { | |
X <- rep(X0, t) | |
for (ti in 2:t) { | |
rd <- randomDisturbance(rd.mu, rd.sigma, stochastic.volatility) | |
X[ti] <- logPriceStep(X[ti - 1], mu, rd, dt) | |
} | |
return(X) | |
} | |
#' @title Generate a discrete price process. | |
#' @description This method exponentiates a log price process and returns it. | |
#' @inheritParams logPriceProcess | |
#' | |
priceProcess <- function(t = 252, X0 = 1.0, mu = 0.1062, | |
rd.mu = 0.0, rd.sigma = 0.45, | |
dt = 0.003968254, stochastic.volatility = TRUE) { | |
return(exp(logPriceProcess(t, log(X0), mu, | |
rd.mu, rd.sigma, | |
dt, stochastic.volatility))) | |
} | |
#' @title Generate an XTS object containing multiple discrete price processes. | |
#' @description This method constructs an xts object which contains multiple | |
#' discrete price processes as simulated using the priceProcess function. These | |
#' price processes are named S1, S2, ... , Sn. | |
#' | |
#' @param n int :: The number of price processes to simulate. | |
#' @inheritParams logPriceProcess | |
#' | |
priceProcesses <- function(n, t = 252, X0 = 1.0, mu = 0.1062, | |
rd.mu = 0.0, rd.sigma = 0.45, | |
dt = 0.003968254, stochastic.volatility = TRUE) { | |
processes <- xtsProcess(priceProcess(t, X0, mu, | |
rd.mu, rd.sigma, | |
dt, stochastic.volatility), "S1") | |
if (n > 1) for (i in 2:n) | |
processes <- merge.xts(processes, | |
xtsProcess(priceProcess(t, X0, mu, rd.mu, rd.sigma, | |
dt, stochastic.volatility), | |
paste("S", i, sep = ''))) | |
return(processes) | |
} | |
#' @title Generate an XTS object containing daily returns for multiple discrete | |
#' price processes. | |
#' | |
#' @description This method constructs an xts object which contains multiple | |
#' daily return processes as simulated using the priceProcess function. These | |
#' price processes are named S1, S2, ... , Sn. | |
#' | |
returnProcesses <- function(n, t = 252, X0 = 1.0, mu = 0.1062, | |
rd.mu = 0.0, rd.sigma = 0.45, dt = 0.003968254, | |
stochastic.volatility = TRUE, method = "discrete") { | |
prices <- priceProcesses(n, t, X0, mu, rd.mu, rd.sigma, dt, stochastic.volatility) | |
return(Return.calculate(prices, method = method)) | |
} | |
#' @title Convert a time series into an xts time series starting from today and | |
#' going forward into the future. | |
#' | |
#' @description This method just converts a time series (a vector of prices or | |
#' returns) into an xts object. | |
#' | |
#' @param ts vector :: The vector containing the time series values. | |
#' @param ts.name character :: The name of the time series. | |
#' | |
xtsProcess <- function(ts, ts.name) { | |
dates <- seq.Date(Sys.Date(), Sys.Date() + (length(ts) - 1), 1) | |
ts.xts <- xts(ts, order.by = dates, col.names = ts.name) | |
colnames(ts.xts) <- ts.name | |
return(ts.xts) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment