Set options and load packages

options(digits=4, width=70)
# install IntroCompFinR package from R-forge
# use install.packages("IntroCompFinR", repos="http://R-Forge.R-project.org")
library(IntroCompFinR)

## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(mvtnorm)
library(PerformanceAnalytics)

## 
## Attaching package: 'PerformanceAnalytics'
## 
## The following object is masked from 'package:graphics':
## 
##     legend

library(zoo)
Sys.setenv(TZ="UTC")

Get data from package IntroCompfinR

# get daily prices on Microsoft from IntroCompFin package
data(msftDailyPrices)
msftPrices = to.monthly(msftDailyPrices, OHLC=FALSE)
# set sample to match book chapter
smpl = "1998-01::2012-05"
msftPrices = msftPrices[smpl]
# calculate returns
msftRetS = Return.calculate(msftPrices, method="simple")
msftRetS = msftRetS[-1]
msftRetC = log(1 + msftRetS)
head(msftRetC, n=3)

##              MSFT
## Feb 1998 0.127864
## Mar 1998 0.054207
## Apr 1998 0.006886

Microsoft cc returns to calibrate Monte Carlo simulation

plot.zoo(msftRetC, main="", 
         lwd=2, col="blue", ylab="Monthly CC Return")
abline(h=0, lwd=2)
abline(h=(mean(msftRetC)+sd(msftRetC)), lty="dashed", lwd=2)
abline(h=(mean(msftRetC)-sd(msftRetC)), lty="dashed", lwd=2)

Mean is 0.0041 and SD is 0.1002

Four panel plot of MSFT cc returns

fourPanelPlot(msftRetC)

Returns have fatter tails than normal and do not show much time dependence.

Simulating observations from the CER model calibrated to MSFT cc returns

mu = 0.004
sd.e = 0.10
n.obs = length(msftRetC)
set.seed(111)
sim.e = rnorm(n.obs, mean=0, sd=sd.e)
sim.ret = mu + sim.e
sim.ret = xts(sim.ret, index(msftRetC))
colnames(sim.ret) = "MSFT.sim"
head(sim.ret, n=3)

##          MSFT.sim
## Feb 1998  0.02752
## Mar 1998 -0.02907
## Apr 1998 -0.02716

Plot simulated returns

plot.zoo(sim.ret, main="", 
         lwd=2, col="blue", ylab="Monthly CC Return")
abline(h=0, lwd=2)
abline(h=(mu+sd.e), lty="dashed", lwd=2)
abline(h=(mu-sd.e), lty="dashed", lwd=2)

Simulated cc returns look like MSFT cc returns except they have constant vol

Four panel plot of simulated returns

fourPanelPlot(sim.ret)

Simulated cc returns share similar stylized facts as MSFT cc returns

Simulated log prices for Msft

sim.p = 2.592 + mu*seq(n.obs) + cumsum(sim.e)
sim.P = exp(sim.p)
# decompose RW into components
sim.sume = zoo(cumsum(sim.e), index(sim.ret))
sim.p = zoo(sim.p, index(sim.ret))
sim.mt = sim.p - sim.sume
sim.P = zoo(sim.P, index(sim.ret))
dataToPlot = merge(sim.sume, sim.mt, sim.p, sim.P)

Plot random walk log prices

plot(dataToPlot, screens=c(1,1,1,2), main="", ylab="",
     lwd=2, col=c("red", "green", "blue", "black"),
     lty=c("dotted", "dashed", "solid", "solid"))
legend(x="center", legend=c("p(t)-E[p(t)]", "E[p(t)]", "log p", "P"),
       lwd=2, col=c("red", "green", "blue", "black"),
       lty=c("dotted", "dashed", "solid", "solid"))

Plot MSFT prices with random walk decomposition

msftLogPrice = log(msftPrices)
p0.msft = as.numeric(log(first(msftPrices)))
mt = p0.msft+ mu*seq(n.obs+1) 
msft.sume = msftLogPrice - mt
dataToPlot = merge(msft.sume, mt, msftLogPrice, msftPrices)
plot.zoo(dataToPlot, screens=c(1,1,1,2), main="", ylab="",
         lwd=2, col=c("red", "green", "blue", "black"),
         lty=c("dotted", "dashed", "solid", "solid"))
legend(x="center", legend=c("p(t)-E[p(t)]", "E[p(t)]", "log p", "P"),
       lwd=2, col=c("red", "green", "blue", "black"),
       lty=c("dotted", "dashed", "solid", "solid"))

Actual log prices look like random walk simulation!

Data for multivariate simulation

# get data from IntroCompFin package
data(sp500DailyPrices, sbuxDailyPrices)
sp500Prices = to.monthly(sp500DailyPrices, OHLC=FALSE)
sbuxPrices = to.monthly(sbuxDailyPrices, OHLC=FALSE)
# set sample to match book chapter
sp500Prices = sp500Prices[smpl]
sbuxPrices = sbuxPrices[smpl]
# calculate returns
sp500RetS = Return.calculate(sp500Prices, method="simple")
sbuxRetS = Return.calculate(sbuxPrices, method="simple")
sp500RetS = sp500RetS[-1]
sbuxRetS = sbuxRetS[-1]
sp500RetC = log(1 + sp500RetS)
sbuxRetC = log(1 + sbuxRetS)
# merged data set
cerRetC = merge(msftRetC, sbuxRetC, sp500RetC)
colnames(cerRetC) = c("MSFT", "SBUX", "SP500")
head(cerRetC, n=3)

##              MSFT    SBUX    SP500
## Feb 1998 0.127864 0.07932 0.068078
## Mar 1998 0.054207 0.13435 0.048738
## Apr 1998 0.006886 0.06096 0.009036

Plot 3 assets

my.panel <- function(...) {
  lines(...)
  abline(h=0)
}
plot.zoo(cerRetC, col="blue", lwd=2, main="", panel=my.panel)

Plot 3 assets

pairs(coredata(cerRetC), col="blue", cex=1.25, pch=16)

Sample descriptive statistics

apply(cerRetC, 2, mean)

##     MSFT     SBUX    SP500 
## 0.004126 0.014657 0.001687

apply(cerRetC, 2, sd)

##    MSFT    SBUX   SP500 
## 0.10021 0.11164 0.04847

cov(cerRetC)

##           MSFT     SBUX    SP500
## MSFT  0.010043 0.003814 0.002998
## SBUX  0.003814 0.012465 0.002476
## SP500 0.002998 0.002476 0.002349

cor(cerRetC)

##         MSFT   SBUX  SP500
## MSFT  1.0000 0.3409 0.6172
## SBUX  0.3409 1.0000 0.4575
## SP500 0.6172 0.4575 1.0000

Multivariate simulation from CER model

# set parameters for simulation
muVec = apply(cerRetC, 2, mean)
covMat = cov(cerRetC)
# multivariate simulation
set.seed(123)
returns.sim = rmvnorm(n.obs, mean=muVec, sigma=covMat)
colnames(returns.sim) = colnames(cerRetC)
returns.sim = zoo(returns.sim, index(cerRetC))
head(returns.sim, n=3)

##              MSFT      SBUX    SP500
## Feb 1998 -0.02303  0.001855  0.05259
## Mar 1998  0.04726  0.054208  0.07676
## Apr 1998  0.01336 -0.125549 -0.03577

Plot simulated returns

plot.zoo(returns.sim, col="blue", lwd=2, main="", panel=my.panel)

Simulated returns look like actual returns except they have constant vol

Plot simulated returns

pairs(coredata(returns.sim), col="blue", cex=1.25, pch=16)

Pairwise correlations in simulated returns match actual returns