Set options and load packages

options(digits=3, 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)

## Warning: package 'mvtnorm' was built under R version 3.2.3

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 data from IntroCompFin package
data(msftDailyPrices, sp500DailyPrices, sbuxDailyPrices)
msftPrices = to.monthly(msftDailyPrices, OHLC=FALSE)
sp500Prices = to.monthly(sp500DailyPrices, OHLC=FALSE)
sbuxPrices = to.monthly(sbuxDailyPrices, OHLC=FALSE)
# set sample to match book chapter
smpl = "1998-01::2012-05"
msftPrices = msftPrices[smpl]
sp500Prices = sp500Prices[smpl]
sbuxPrices = sbuxPrices[smpl]
# calculate returns
msftRetS = na.omit(Return.calculate(msftPrices, method="simple"))
sp500RetS = na.omit(Return.calculate(sp500Prices, method="simple"))
sbuxRetS = na.omit(Return.calculate(sbuxPrices, method="simple"))
msftRetC = log(1 + msftRetS)
sp500RetC = log(1 + sp500RetS)
sbuxRetC = log(1 + sbuxRetS)
# merged data set
cerRetS = merge(msftRetS, sbuxRetS, sp500RetS)
cerRetC = merge(msftRetC, sbuxRetC, sp500RetC)
colnames(cerRetS) = colnames(cerRetC) = c("MSFT", "SBUX", "SP500")
head(cerRetC, n=3)

##             MSFT   SBUX   SP500
## Feb 1998 0.12786 0.0793 0.06808
## Mar 1998 0.05421 0.1343 0.04874
## Apr 1998 0.00689 0.0610 0.00904

Plot 3 assets

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

Plot 3 assets

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

CER model parameter estimates = sample descriptive statistics

# Estimate mean (mu)
muhat = apply(cerRetC,2,mean)
muhat

##    MSFT    SBUX   SP500 
## 0.00413 0.01466 0.00169

# Estimate variance (sigma^2)
sigma2hat = apply(cerRetC,2,var)
sigma2hat

##    MSFT    SBUX   SP500 
## 0.01004 0.01246 0.00235

# Estimate volatility (sigma)
sigmahat = apply(cerRetC,2,sd)
sigmahat

##   MSFT   SBUX  SP500 
## 0.1002 0.1116 0.0485

Show mean-volatility tradeoff

plot(sigmahat, muhat, col="blue", pch=16, cex=2,
     ylim=c(0, 0.015), xlim=c(0, 0.12))
text(sigmahat, muhat, labels=colnames(cerRetC), pos=4)

Assets with highest means also have highest volatilities

Estimate covariance and correlation matrices

covmat = var(cerRetC)
covmat

##          MSFT    SBUX   SP500
## MSFT  0.01004 0.00381 0.00300
## SBUX  0.00381 0.01246 0.00248
## SP500 0.00300 0.00248 0.00235

cormat = cor(cerRetC)
cormat

##        MSFT  SBUX SP500
## MSFT  1.000 0.341 0.617
## SBUX  0.341 1.000 0.457
## SP500 0.617 0.457 1.000

Extract unique values from covariance and correlation matrices

covhat = covmat[lower.tri(covmat)]
rhohat = cormat[lower.tri(cormat)]
names(covhat) <- names(rhohat) <- 
c("msft,sbux","msft,sp500","sbux,sp500")
covhat

##  msft,sbux msft,sp500 sbux,sp500 
##    0.00381    0.00300    0.00248

rhohat

##  msft,sbux msft,sp500 sbux,sp500 
##      0.341      0.617      0.457

Estimated standard errors for means

n.obs = nrow(cerRetC)
seMuhat = sigmahat/sqrt(n.obs)
cbind(muhat, seMuhat)

##         muhat seMuhat
## MSFT  0.00413 0.00764
## SBUX  0.01466 0.00851
## SP500 0.00169 0.00370

Estimated SE values are large compared to estimates. Do not estimate mean values very well.

Estimated standard errors for variances and standard deviations

seSigma2hat = sigma2hat/sqrt(n.obs/2)
seSigmahat = sigmahat/sqrt(2*n.obs)
cbind(sigma2hat, seSigma2hat, sigmahat, seSigmahat)

##       sigma2hat seSigma2hat sigmahat seSigmahat
## MSFT    0.01004    0.001083   0.1002    0.00540
## SBUX    0.01246    0.001344   0.1116    0.00602
## SP500   0.00235    0.000253   0.0485    0.00261

Estimated SE values for volatility are small compared to estimates. Estimate volatility well.

Estimated standard errors for correlations

seRhohat = (1-rhohat^2)/sqrt(n.obs)
cbind(rhohat, seRhohat)

##            rhohat seRhohat
## msft,sbux   0.341   0.0674
## msft,sp500  0.617   0.0472
## sbux,sp500  0.457   0.0603

Estimated standard errors for correlations are somewhat small compared to estimates. Notice how standard error is smallest for the largest correlation estimate.

Approximate 95% confidence intervals for \(\mu\)

lowerMu = muhat - 2*seMuhat
upperMu = muhat + 2*seMuhat
widthMu = upperMu - lowerMu
cbind(lowerMu, upperMu, widthMu)

##        lowerMu upperMu widthMu
## MSFT  -0.01116 0.01941  0.0306
## SBUX  -0.00237 0.03168  0.0341
## SP500 -0.00570 0.00908  0.0148

Wide 95% confidence intervals for the mean implies imprecise estimates. Notice that the confidence interval width is smallest for the S&P 500 index (why?)

Approximate 95% confidence intervals for \(\sigma^2\) and \(\sigma\)

# confidence intervals for sigma^2
lowerSigma2 = sigma2hat - 2*seSigma2hat
upperSigma2 = sigma2hat + 2*seSigma2hat
widthSigma2 = upperSigma2 - lowerSigma2
cbind(lowerSigma2, upperSigma2, widthSigma2)

##       lowerSigma2 upperSigma2 widthSigma2
## MSFT      0.00788     0.01221     0.00433
## SBUX      0.00978     0.01515     0.00538
## SP500     0.00184     0.00286     0.00101

# confidence intervals for sigma
lowerSigma = sigmahat - 2*seSigmahat
upperSigma = sigmahat + 2*seSigmahat
widthSigma = upperSigma - lowerSigma
cbind(lowerSigma, upperSigma, widthSigma)

##       lowerSigma upperSigma widthSigma
## MSFT      0.0894     0.1110     0.0216
## SBUX      0.0996     0.1237     0.0241
## SP500     0.0432     0.0537     0.0105

Confidence intervals for sigma are narrow. Estimates are precise.

Approximate 95% confidence intervals for \(\rho\)

lowerRho = rhohat - 2*seRhohat
upperRho = rhohat + 2*seRhohat
widthRho = upperRho - lowerRho
cbind(lowerRho, upperRho, widthRho)

##            lowerRho upperRho widthRho
## msft,sbux     0.206    0.476    0.270
## msft,sp500    0.523    0.712    0.189
## sbux,sp500    0.337    0.578    0.241

Confidence intervals are not too wide and contain all positive values. Hence, estimates are moderately precise.

Stylized facts for the estimation of CER model parameters

The expected return is not estimated very precisely
- Large standard errors relative to size of mean estimates
Standard deviations and correlations are estimated more precisely than the expected return

Monte Carlo Simulation Loop

mu = 0.05
sigma = 0.10
n.obs = 100
n.sim = 1000
set.seed(111)
sim.means = rep(0,n.sim)  # initialize vectors
mu.lower = rep(0,n.sim)
mu.upper = rep(0,n.sim)
qt.975 = qt(0.975, n.obs-1)
for (sim in 1:n.sim) {
    sim.ret = rnorm(n.obs,mean=mu,sd=sigma)
    sim.means[sim] = mean(sim.ret)
    se.muhat = sd(sim.ret)/sqrt(n.obs)
    mu.lower[sim] = sim.means[sim]-qt.975*se.muhat
    mu.upper[sim] = sim.means[sim]+qt.975*se.muhat
}

First 10 simulated samples from CER model

Histogram of 1000 Monte Carlo estimates of \(\mu\)

hist(sim.means, col="cornflowerblue", ylim=c(0,40), 
     main="", xlab="muhat", probability=TRUE)
abline(v=mean(sim.means), col="white", lwd=4, lty=2)
# overlay normal curve
x.vals = seq(0.02, 0.08, length=100)
lines(x.vals, dnorm(x.vals, mean=mu, sd=sigma/sqrt(100)), col="orange", lwd=2)

Approximate expected value, bias, SE, and 95% CI coverage probability

mean(sim.means)       # MC expected value

## [1] 0.0497

mean(sim.means) - mu  # MC bias

## [1] -0.00031

sd(sim.means)         # MC SE

## [1] 0.0104

sigma/sqrt(n.obs)     # true SE(muhat)

## [1] 0.01

# 95% CI coverage probability
in.interval = mu >= mu.lower & mu <= mu.upper
sum(in.interval)/n.sim

## [1] 0.931

Histograms of 1000 Monte Carlo estimates of \(\sigma^2\) and \(\sigma\)

Approximate expected value, bias, SE, and 95% CI coverage probability

# Evaluation of bias
mean(sim.sigma2)           # true sigma^2 = 0.01

## [1] 0.00999

mean(sim.sigma2) - sigma^2 # MC bias for sigma^2

## [1] -9.87e-06

mean(sim.sigma)            # true sigma = 0.1

## [1] 0.0997

mean(sim.sigma) - sigma    # MC bias for sigma

## [1] -0.000278

# Evaluation of SE
sd(sim.sigma2)         # MC SE for sigma^2

## [1] 0.00135

sigma^2/sqrt(n.obs/2)  # Asymptotic SE

## [1] 0.00141

sd(sim.sigma)          # MC SE for sigma

## [1] 0.00676

sigma/sqrt(2*n.obs)    # Asymptotic SE

## [1] 0.00707

# Evaluation of 95% Confidence Interval for sigma2
in.interval = sigma^2 >= sigma2.lower & sigma^2 <= sigma2.upper
sum(in.interval)/n.sim

## [1] 0.951

# Evaluation of 95% Confidence Interval for sigma
in.interval = sigma >= sigma.lower & sigma <= sigma.upper
sum(in.interval)/n.sim

## [1] 0.963

Monte Carlo SE values are close to asymptotic SE values.
Monte Carlo SE values are more accurate (provided number of simulations is large)
Monte Carlo 95% confidence interval coverage is close to 95%

Histogram of 1000 Monte Carlo estimates of \(\rho\)

True value is 0.75

Approximate expected value, bias, SE, and 95% CI coverage probability

mean(sim.corrs)          # true value is 0.75

## [1] 0.747

mean(sim.corrs) - rho12  # bias

## [1] -0.00315

sd(sim.corrs)            # MC SE

## [1] 0.045

(1-rho12^2)/sqrt(n.obs)  # Asymptotic SE

## [1] 0.0437

# Evaluation of 95% Confidence Interval for rho
in.interval = (rho12 >= rho.lower) & (rho12 <= rho.upper)
sum(in.interval)/n.sim

## [1] 0.952

Estimating simple return quantiles from the CER Model

n.obs = length(msftRetC)
muhatS = colMeans(cerRetS)
sigmahatS = apply(cerRetS, 2, sd)
qhat.05 = muhatS + sigmahatS*qnorm(0.05)
qhat.01 = muhatS + sigmahatS*qnorm(0.01)
seQhat.05 = (sigmahatS/sqrt(n.obs))*sqrt(1 + 0.5*qnorm(0.05)^2)
seQhat.01 = (sigmahatS/sqrt(n.obs))*sqrt(1 + 0.5*qnorm(0.01)^2)
seQhat.001 = (sigmahatS/sqrt(n.obs))*sqrt(1 + 0.5*qnorm(0.001)^2)
qhat.001 = muhatS + sigmahatS*qnorm(0.001)
cbind(qhat.05, seQhat.05, qhat.01, seQhat.01, qhat.001, seQhat.001)

##       qhat.05 seQhat.05 qhat.01 seQhat.01 qhat.001 seQhat.001
## MSFT  -0.1578   0.01187  -0.227   0.01490   -0.304    0.01860
## SBUX  -0.1587   0.01276  -0.233   0.01602   -0.316    0.02000
## SP500 -0.0758   0.00559  -0.108   0.00702   -0.145    0.00876

SE values for 1% quantiles and bigger than SE values for 5% quantiles

95% Confidence Intervals for simple return quantiles

# 5% quantiles
lowerQhat.05 = qhat.05 - 2*seQhat.05
upperQhat.05 = qhat.05 + 2*seQhat.05
widthQhat.05 = upperQhat.05 - lowerQhat.05
cbind(lowerQhat.05, upperQhat.05, widthQhat.05)

##       lowerQhat.05 upperQhat.05 widthQhat.05
## MSFT       -0.1815      -0.1341       0.0475
## SBUX       -0.1842      -0.1331       0.0511
## SP500      -0.0869      -0.0646       0.0224

# 1% quantiles
lowerQhat.01 = qhat.01 - 2*seQhat.01
upperQhat.01 = qhat.01 + 2*seQhat.01
widthQhat.01 = upperQhat.01 - lowerQhat.01
cbind(lowerQhat.01, upperQhat.01, widthQhat.01)

##       lowerQhat.01 upperQhat.01 widthQhat.01
## MSFT        -0.257      -0.1972       0.0596
## SBUX        -0.265      -0.2010       0.0641
## SP500       -0.122      -0.0943       0.0281

# 95% CI for .1% quantile
lowerQhat.001 = qhat.001 - 2*seQhat.001
upperQhat.001 = qhat.001 + 2*seQhat.001
widthQhat.001 = upperQhat.001 - lowerQhat.001
cbind(lowerQhat.001, upperQhat.001, widthQhat.001)

##       lowerQhat.001 upperQhat.001 widthQhat.001
## MSFT         -0.342        -0.267        0.0744
## SBUX         -0.356        -0.276        0.0800
## SP500        -0.162        -0.127        0.0350

95% confidence intervals are somewhat large

Estimation of CER Model

Set options and load packages

Get data from package IntroCompfinR

Plot 3 assets

Plot 3 assets

CER model parameter estimates = sample descriptive statistics

Show mean-volatility tradeoff

Estimate covariance and correlation matrices

Extract unique values from covariance and correlation matrices

Estimated standard errors for means

Estimated standard errors for variances and standard deviations

Estimated standard errors for correlations

Approximate 95% confidence intervals for \(\mu\)

Approximate 95% confidence intervals for \(\sigma^2\) and \(\sigma\)

Approximate 95% confidence intervals for \(\rho\)

Stylized facts for the estimation of CER model parameters

Monte Carlo Simulation Loop

First 10 simulated samples from CER model

Histogram of 1000 Monte Carlo estimates of \(\mu\)

Histograms of 1000 Monte Carlo estimates of \(\sigma^2\) and \(\sigma\)

Approximate expected value, bias, SE, and 95% CI coverage probability

Histogram of 1000 Monte Carlo estimates of \(\rho\)

Approximate expected value, bias, SE, and 95% CI coverage probability

Estimating simple return quantiles from the CER Model

95% Confidence Intervals for simple return quantiles

Histograms of 1000 Monte Carlo estimates of \(\hat{q}_{\alpha}\) and \(VaR_{\alpha}\)