R/stock_price_simulation.R

In mrepsilon/PawelKawskiPackage: Stock Analytics Package

#' Generate stock prices
#'
#' This function plots a chosen number of random trajectories for stock marked prices of chosen asset.
#' Function on its own is providing data from google finance. Stochastic process which is used in function
#' to modeling future prices is a Black Scholes model. Simulation values are aproximations of process: \cr
#' \eqn{dX_t = a * X_t dt + sigma * X_t dW_t} , \cr
#' where W is a standard Wiener process.
#'
#' @param symbol Symbol of a stock exchange data.
#' @param days Number of days for which simulation is to be performed.
#' @param simulations Number of simulated trajectories which has to be done.
#' @param price_type Type of price : "close" or "open".
#' @param driftless Logical value which can set dift on 0 in model. Default value is FALSE.
#' @return A plot of simulated process trajectories and list with 2 frames ( [[1]] - historical data and [[2]] - simulated trajectories).
#' @examples
#' stock_price_simulation("GOOG",100,200,"open",TRUE)
#' stock_price_simulation("YHOO",100,200)
#' stock_price_simulation("GOOG",100,50,"close")
#' @importFrom graphics lines plot
#' @importFrom stats na.omit quantile rnorm sd
#' @export

stock_price_simulation <- function(symbol, days, simulations, price_type = "close", driftless = FALSE){

  assertthat::assert_that(is.character(symbol))
  assertthat::assert_that(days%%1==0 & days>=1, msg="Number of days must be integer value and greater then 0")
  assertthat::assert_that(price_type=="close" | price_type=="open")
  assertthat::assert_that(is.logical(driftless), msg="Do you want simulate process of prices with drift or not? It's logical value.")
  assertthat::assert_that(simulations%%1==0, msg="Number of simulations must be integer value and greater then 0")

  # basic parameters
  t <- days/365
  dt <- 1/365
  seq1 <- seq(0, t, by = dt)
  N <- length(seq1)
  m <- simulations

  # get data from google finance
  options(warn = -1)
  data <- quantmod::getSymbols (symbol, src = "google", env = NULL)
  options(warn = 1)

  # test that the columns are in set order
  assertthat::assert_that(grepl("Open",colnames(data)[1]))
  assertthat::assert_that(grepl("Close",colnames(data)[4]))

  if (price_type=="close"){data_col <- 4}
  if (price_type=="open"){data_col <- 1}

  # simple parameters estimation from data
  differences <- diff(data[, data_col])

  drift <- if(driftless==F) {
                            mean(differences, na.rm = T)
           }else{drift <- 0}

  risk <- 2 # it is only assumption used because this model is builded under real world measure and Black Scholes is constructed under (martingale) risk neutral measure
  volatility <- sd(as.numeric(differences), na.rm = T)/mean(data[,data_col],na.rm = T) * risk


  # simulation using Black Scholes model for stock prices.
  Y <- matrix(rep(0, N * m), ncol = N, nrow = m)
  dB.matrix <- matrix(rep(0, (N - 1) * m), ncol = (N - 1), nrow = m)

  for(j in 1:m){
    X <- c()
    X[1] <- as.numeric(data[nrow(data), data_col])
    dB <- sqrt(dt) * rnorm(N - 1, mean = 0, sd = 1)
    dB.matrix[j, ] <- dB
    for(i in 1:(N - 1)){
      X[i + 1] <- X[i] + X[i] * drift * dt + X[i] * volatility * dB[i]
    }
    Y[j, ] <- X
  }

  #plot
  darkcols <- RColorBrewer::brewer.pal(8, "Dark2")
  if(m>8){darkcols<-rep(darkcols,floor(m/8)+1)}

  ##create sequence of dates (it is not important that they are not business days for our plot)
  days_seq <- seq(zoo::index(data)[nrow(data)], (zoo::index(data)[nrow(data)] + days) , by="day" )

  #filter data to about last year (more visible plot)
  truncate_date <- nrow(data)-300
  if(nrow(data)>301){data <- data[truncate_date : nrow(data), ]}

  simulated_trajectories <- Y

  #plot data and forecasts
 plot(zoo::index(data), as.numeric(data[,data_col]),
        xlim = c(zoo::index(data)[1], zoo::index(data)[nrow(data)] + days),
        ylim=c(min(min(as.numeric(data[,data_col])),min(Y)), max(max(as.numeric(data[,data_col])), max(Y) ) ),
        ylab = "Price", xlab = "Date",
        type="l", main = "Simulated trajectories")
  for(j in 2:m){
    lines(days_seq, Y[j,], col = darkcols[j])
  }
   return(list(data,simulated_trajectories))
}