#######################
## Only for educational purposes, not intended as an investment tool.
## Covers the base case, but can be adapted to more assets or different
## deliveries or time steps.
## This code has been possible thanks to the generosity of the:
## http://stackoverflow.com/ community
## Thanks to http://www.quantmod.com/ excellent package
## Contact:
## https://www.researchgate.net/profile/Mariano_Mendez
#######################
##### Just in case...
# install.packages("quantmod")
# install.packages("ggthemes")
# install.packages("scales")
# install.packages("MASS")
# install.packages("ggplot2")
# install.packages("reshape2")
# install.packages("corrplot")
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: TTR
## Version 0.4-0 included new data defaults. See ?getSymbols.
library(ggthemes)
library(scales)
library(MASS)
library(ggplot2)
library(reshape2)
library(corrplot)
rm(list = ls())
set.seed(233)
### Product especific information
timeSteps <- 1 # number of time steps inside periods
periods <- 10 * timeSteps # number of periods in the simulation
periodsPerYear <- 1 # one year
sigmaTime <- 52 # Weeks semesters, adjust sigma to time of simulation
t <- periods + 1 # number of periods plus 1
dt <- 1 # time interval 1
paths <- 1e4
pathsGraph <- 190 # Only small number of paths for graph simulation
### Information for retrieving asset data
ticker <- c("JPM","C") ### finance.yahoo.com tickers
assets <- length(ticker)
from <- '2004-01-01'
to <- '2008-05-22' # Historical data: yyyy-mm-dd
rf <- c(0.0388, 0.0388)/(periodsPerYear) # rf adjusted by periods per year
strike <- 40
# Labels for the simulation graph
ValuationDate <- rep(to,t)
ValuationDate <- as.POSIXlt(to)
ValuationDate$year <-ValuationDate$year+c(0:periods)
ValuationDate <- format(ValuationDate,"%Y")
### Obtaining prices, all at de same time
getSymbols(ticker,from=from,to=to)
## As of 0.4-0, 'getSymbols' uses env=parent.frame() and
## auto.assign=TRUE by default.
##
## This behavior will be phased out in 0.5-0 when the call will
## default to use auto.assign=FALSE. getOption("getSymbols.env") and
## getOptions("getSymbols.auto.assign") are now checked for alternate defaults
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for more details.
## [1] "JPM" "C"
ClosePrices <- do.call(cbind,lapply(ticker,function(x) Cl(get(x))))
ClosePrices <- na.locf(ClosePrices) # get prices but getting rid of NA
AdjustedPrices <- do.call(cbind,lapply(ticker,function(x) Ad(get(x))))
AdjustedPrices <- na.locf(AdjustedPrices) # get prices but getting rid of NA
ReturnPrices <- do.call(cbind,lapply(ticker,function(x) weeklyReturn(get(x))))
ReturnPrices <- na.locf(ReturnPrices) # get prices but getting rid of NA
colnames(ReturnPrices) <- ticker
colnames(AdjustedPrices) <- ticker
### Calculating dividend yield, because of quantmod we need to ge one at a time
### Because Citibank Reverse Split we use adjusted close
Dividends <- list()
CheckDividends <- list()
for(i in seq_along(ticker)){
AdcloseP <- get(ticker[[i]])[,6] #get adjusted closing prices
div <- getDividends(ticker[i], from = from, to = to, auto.assign = FALSE)
if(nrow(div!=0)){
temp <- merge.xts(AdcloseP,div,join='right')
temp$yld <- temp[,2]/temp[,1]
temp1 <- apply.yearly(temp$yld,sum)
temp1 <- last(temp1, '-1 year')
CheckDividends[[i]] <- temp1
temp1 <- colMeans(temp1)}
else{temp1 <- 0}
Dividends <- c(Dividends,temp1)
}
## some cleaning...
AdcloseP <- NULL; div <- NULL; temp <- NULL; temp1 <- NULL
### Gathering data for the simulation of GBM
Dividends <- as.numeric(Dividends)
GraphDividends <-data.frame(ticker, Dividends) # Yearly div for graph
Dividends <- Dividends/(periodsPerYear)
## Choose model with or without dividends
#Dividends <- 0
#### historical volatility estimation
sigma <- apply(ReturnPrices, 2, sd)
sigma <- sigma*sqrt(sigmaTime) # from weekly to semester sigma
correlation <- cor(ReturnPrices)
# Setting initial prices for the simulation
InitialPrices <- ClosePrices[nrow(ClosePrices),]
colnames(InitialPrices) <- ticker
S0 <- ClosePrices[nrow(ClosePrices),]
S0 <- matrix(rep(S0,each=paths),nrow=paths)
# First and second terms of GBM
g1 <- matrix(rep(rf-Dividends-((sigma^2)/2)*dt,each=paths),nrow = paths)
g2 <- matrix(rep(sigma*sqrt(dt),each=paths),nrow = paths)
###############
#### Uncorrelated noise matrix
wiener <- rnorm(paths*t*assets)
dim(wiener) <- c(paths,t,assets)
#################
# ### Correlated noise matrix
# wiener <- mvrnorm(paths*t,mu=rep(0, assets),Sigma=correlation,empirical=T)
# dim(wiener) <- c(paths, t, assets)
####################
# Preasigment of variables for Loops
St <- array(0,c(paths,t,assets))
returns <- list()
stock <- list()
# Generating simulated prices
for(i in 1:assets){
St[,,i] <- exp(g1[,i]+wiener[,,i]*g2[,i])
St[,1,] <- S0
St[,,i] <- t(apply(St[,,i],1,cumprod))
stock[[i]] <- St[,,i]
}
## Call valuation by simulation
last_price <- c(t,(t+t))
valuation <- data.frame(stock)
valuation <- valuation[last_price]
call_sim <- valuation-strike
call_sim[call_sim<0] <- 0
call_sim <- exp(-rf*periods)*colMeans(call_sim)
## BSM Valuation
InitialPrices <- as.numeric(InitialPrices)
d1 <- (log(InitialPrices/strike)+(rf-Dividends+sigma^2/2)*periods)/(sigma*sqrt(periods))
d2 <- d1 - sigma*sqrt(periods)
call_BSM <- InitialPrices*exp(-Dividends*periods)*pnorm(d1)-strike*exp(-rf*periods)*pnorm(d2)
### Binding Output
call_Prices <- rbind(InitialPrices, strike, sigma, rf,
Dividends, call_sim, call_BSM)
colnames(call_Prices) <- ticker
##################
### GRAPHS########
##################
## Historical Prices Graph closing prices
df<-data.frame(Date=index(ClosePrices),coredata(ClosePrices))
df_melt = melt(df, id.vars = 'Date')
ggplot(df_melt, aes(x = Date, y = value)) +
geom_line() +
facet_wrap(~ variable, scales = 'free_y', ncol = 1) +
ylab("Closing Prices") +
theme(axis.title.x=element_blank())
## Dividends Graph
ggplot(GraphDividends,aes(x=ticker,y=Dividends,fill=factor(ticker))) +
geom_point(size=3) +
ggtitle("Annual dividend rate") +
scale_y_continuous(breaks = scales::pretty_breaks(n = 15),labels=percent) +
theme_economist() +
scale_colour_economist() +
theme(legend.position="none",axis.title.x = element_blank(),axis.title.y=element_blank())
###Correlation
corrplot(correlation, addCoef.col = "black", type="upper")
#
## Simulation
# Row binding the different stocks for the "pathsGraph" number of simulated paths
# do.call calls the function
stock <- do.call( rbind, lapply(stock,'[',1:pathsGraph,))
group <- rep(ticker,each=pathsGraph)
stock <- data.frame(id=seq_along(group),group,stock)
stockm <- melt(stock,id.vars = c("id","group"))
ggplot(stockm, aes(variable, value, group = id)) +
geom_line(aes(colour = factor(id)),alpha = 0.5) +
facet_wrap(~ group, scales="free_y", ncol = 1) +
theme(legend.position="none") +
scale_x_discrete(labels=as.character.Date(ValuationDate),expand = c(0,0)) +
xlab("Time") + ylab(paste0("Sample of ",pathsGraph," paths for each asset"))+
geom_hline(aes(yintercept=strike),color="black",size=0.8, linetype="dashed")
call_Prices
## JPM C
## InitialPrices 43.04999900 21.72000000
## strike 40.00000000 40.00000000
## sigma 0.25336406 0.26496261
## rf 0.03880000 0.03880000
## Dividends 0.04483412 0.04602221
## call_sim 8.92754864 1.85333964
## call_BSM 8.68370747 1.79619654