# econ424lab5.r			script file for econ 424 lab4 calculations
#
# author: Eric Zivot
# created: October 20, 2003
# revised: November 1, 2009
#
# comments:
# Data for the lab are in the csv file econ424lab4returns.csv, which contains 
# monthly continuously compounded returns on Vanguard long term bond index fund 
# (VBLTX), Fidelity Magellan stock mutual fund (FMAGX), and General Motors 
# stock (GM) 
#
options(digits=4)
options(chmhelp=TRUE)

library(PerformanceAnalytics)
library("zoo")
library("mvtnorm")
library("boot")

# load data
# it is assumed that the file is in the location 
# C:/Users/ezivot/Documents/classes/econ424/fall2009
# change this location for your application.

lab4.df = read.csv("C:/Users/ezivot/Documents/classes/econ424/fall2009/econ424lab4returns.csv",
                  stringsAsFactors=F)
colnames(lab4.df)
head(lab4.df)

#
# 2. Create zoo object from data and dates in lab4.df. This object is only used
# for the time series plots. See the document "Working with Time Series in R" on the
# class webpage for more details on zoo objects
#    

lab4.z = zoo(x=lab4.df[, -1], 
             order.by=as.yearmon(lab4.df[, 1], format="%b-%y"))
start(lab4.z)
end(lab4.z)
colnames(lab4.z) 
head(lab4.z)


#
# 3. Create timePlots of data
#

plot(lab4.z,col="blue", lwd=2)
plot(lab4.z, plot.type="single", lty=1:3, col=1:3, lwd=2)
legend(x="bottomleft", legend=colnames(lab4.z), lty=1:3, col=1:3, lwd=2)
abline(h=0)
title("Monthly cc returns")


#
# 4. Create matrix of return data and compute pairwise scatterplots
#

ret.mat = coredata(lab4.z)
colnames(ret.mat)
head(ret.mat)
VBLTX = ret.mat[,"VBLTX"]
FMAGX = ret.mat[,"FMAGX"]
GM = ret.mat[,"GM"]
pairs(ret.mat, col="blue")

#
# 5. Compute estimates of CER model parameters
#
muhat.vals = apply(ret.mat, 2, mean)
muhat.vals
sigma2hat.vals = apply(ret.mat, 2, var)
sigma2hat.vals
sigmahat.vals = apply(ret.mat, 2, sd)
sigmahat.vals
cov.mat = var(ret.mat)
cov.mat
cor.mat = cor(ret.mat)
cor.mat
covhat.vals = cov.mat[lower.tri(cov.mat)]
rhohat.vals = cor.mat[lower.tri(cor.mat)]
names(covhat.vals) <- names(rhohat.vals) <- 
c("VBLTX,FMAGX","VBLTX,GM","FMAGX,GM")
covhat.vals
rhohat.vals

# summarize the CER model estimates
cbind(muhat.vals,sigma2hat.vals,sigmahat.vals)
cbind(covhat.vals,rhohat.vals)

# plot mean vs. sd values
plot(sigmahat.vals, muhat.vals, pch=1:3, cex=2, col=1:3, 
     ylab = "mean", xlab="sd (risk)", ylim=c(0, 0.01))
legend(x="topleft", legend=names(muhat.vals), pch=1:3, col=1:3, cex=1.5)     

#
# 6. Compute stndard errors for estimated parameters
#

# compute estimated standard error for mean
nobs = nrow(ret.mat)
nobs
se.muhat = sigmahat.vals/sqrt(nobs)
se.muhat
# show estimates with SE values underneath
rbind(muhat.vals,se.muhat)

# compute approx 95% confidence intervals
mu.lower = muhat.vals - 2*se.muhat
mu.upper = muhat.vals + 2*se.muhat
cbind(mu.lower,mu.upper)

# compute estimated standard errors for variance and sd
se.sigma2hat = sigma2hat.vals/sqrt(nobs/2)
se.sigma2hat
se.sigmahat = sigmahat.vals/sqrt(2*nobs)
se.sigmahat

rbind(sigma2hat.vals,se.sigma2hat)
rbind(sigmahat.vals,se.sigmahat)

# compute approx 95% confidence intervals
sigma2.lower = sigma2hat.vals - 2*se.sigma2hat
sigma2.upper = sigma2hat.vals + 2*se.sigma2hat
cbind(sigma2.lower,sigma2.upper)

sigma.lower = sigmahat.vals - 2*se.sigmahat
sigma.upper = sigmahat.vals + 2*se.sigmahat
cbind(sigma.lower,sigma.upper)

# compute estimated standard errors for correlation
se.rhohat = (1-rhohat.vals^2)/sqrt(nobs)
se.rhohat
rbind(rhohat.vals,se.rhohat)

# compute approx 95% confidence intervals
rho.lower = rhohat.vals - 2*se.rhohat
rho.upper = rhohat.vals + 2*se.rhohat
cbind(rho.lower,rho.upper)


#
# 7. Compute 5% and 1% Value at Risk
#

# function to compute Value-at-Risk
# note: default values are selected for 
# the probability level (p) and the initial
# wealth (w). These values can be changed
# when calling the function. Highlight the entire
# function, right click and select run line or selection
Value.at.Risk = function(x,p=0.05,w=100000) {
	x = as.matrix(x)
	q = apply(x, 2, mean) + apply(x, 2, sd)*qnorm(p)
	VaR = (exp(q) - 1)*w
	VaR
}

# 5% and 1% VaR estimates based on W0 = 100000

Value.at.Risk(ret.mat,p=0.05,w=100000)
Value.at.Risk(ret.mat,p=0.01,w=100000)


#
# 8. Evaluate bias and SE formulas using Monte Carlo
#

# generate 1000 samples from CER and compute sample statistics

mu = muhat.vals["FMAGX"]
sd = sigmahat.vals["FMAGX"]
n.obs = 60
set.seed(123)
n.sim = 1000

sim.means = rep(0,n.sim)
sim.vars = rep(0,n.sim)
sim.sds = rep(0,n.sim)
for (sim in 1:n.sim) {
	sim.ret = rnorm(n.obs,mean=mu,sd=sd)
	sim.means[sim] = mean(sim.ret)
	sim.vars[sim] = var(sim.ret)
	sim.sds[sim] = sqrt(sim.vars[sim])
}

par(mfrow=c(2,2))
hist(sim.means,xlab="mu hat", col="slateblue1")
abline(v=mu, col="white", lwd=2)
hist(sim.vars,xlab="sigma2 hat", col="slateblue1")
abline(v=sd^2, col="white", lwd=2)
hist(sim.sds,xlab="sigma hat", col="slateblue1")
abline(v=sd, col="white", lwd=2)
par(mfrow=c(1,1))

# 
# 9. compute MC estimates of bias and SE
#

c(mu, mean(sim.means))
mean(sim.means) - mu
c(sd^2, mean(sim.vars))
mean(sim.vars) - sd^2
c(sd, mean(sim.sds))
mean(sim.sds) - sd

# compute MC SE value and compare to SE calculated from simulated data

c(se.muhat["FMAGX"], sd(sim.means))
c(se.sigma2hat["FMAGX"], sd(sim.vars)) 
c(se.sigmahat["FMAGX"], sd(sim.sds))

# 
# 10. simulate data for three asset returns
#

mu = muhat.vals
cov = cov.mat
nobs = nrow(ret.mat)
set.seed(123)
sim.e = rmvnorm(nobs, mean=rep(0,3), sigma=cov.mat)
sim.ret = mu + sim.e
colnames(sim.ret) = paste(colnames(ret.mat),".sim",sep="")

# plot data and scatterplots
ts.plot(sim.ret,main="Simulated return data", lty=1:3, col=1:3, lwd=2)
abline(h=0)
legend(x="topleft",legend=colnames(sim.ret), lty=1:3, col=1:3, lwd=2)
pairs(sim.ret, col="blue")

# compute covariances and correlations
cov.mat.sim = var(sim.ret)
cor.mat.sim = cor(sim.ret)
covhat.sim.vals = cov.mat.sim[lower.tri(cov.mat)]
rhohat.sim.vals = cor.mat.sim[lower.tri(cor.mat)]
names(covhat.sim.vals) <- names(rhohat.sim.vals) <- 
c("VBLTX.sim,FMAGX.sim","VBLTX.sim,GM.sim","FMAGX.sim,GM.sim")
covhat.sim.vals
rhohat.sim.vals
# compute estimated standard errors for correlation
se.rhohat.sim = (1-rhohat.sim.vals^2)/sqrt(nobs)
se.rhohat.sim
# show estimates with se values underneath
rbind(rhohat.sim.vals, se.rhohat.sim)

# compute approx 95% confidence intervals
rho.sim.lower = rhohat.sim.vals - 2*se.rhohat.sim
rho.sim.upper = rhohat.sim.vals + 2*se.rhohat.sim
cbind(rho.sim.lower,rho.sim.upper)

# generate 1000 samples from CER and compute correlations
n.obs = 60
n.sim = 1000
set.seed(123)
sim.corrs = matrix(0,n.sim,3)	# initialize vectors
colnames(sim.corrs) = c("vbltx,fmagx","vbltx,gm","fmagx,gm")

for (sim in 1:n.sim) {
	sim.ret = rmvnorm(n.obs,mean=muhat.vals, sigma=cov.mat)
	cor.mat = cor(sim.ret)
	sim.corrs[sim,] = cor.mat[lower.tri(cor.mat)]
}

par(mfrow=c(2,2))
	hist(sim.corrs[,1], xlab="vbltx,fmagx", col="slateblue1", main="")
	abline(v=rhohat.vals[1], col="white", lwd=2)
	hist(sim.corrs[,2], xlab="vbltx,gm", col="slateblue1", main="")
		abline(v=rhohat.vals[2], col="white", lwd=2)
	hist(sim.corrs[,3], xlab="fmagx,gm", col="slateblue1", main="")
		abline(v=rhohat.vals[3], col="white", lwd=2)
par(mfrow=c(1,1))

# MC means and standard deviations
apply(sim.corrs, 2, mean)
# bias
apply(sim.corrs, 2, mean) - rhohat.vals
apply(sim.corrs, 2, sd)
se.rhohat

#
# 11. bootstrapping SE for mean, variance, sd and correlation
#

?boot
# note: boot requires user-supplied functions that take
# two arguments: data and an index. The index is created
# by the boot function and represents random resampling with
# replacement

# function for bootstrapping sample mean
mean.boot = function(x, idx) {
# arguments:
# x 		data to be resampled
# idx		vector of scrambled indices created by boot() function
# value:
# ans		mean value computed using resampled data
     ans = mean(x[idx])
     ans
}

VBLTX.mean.boot = boot(VBLTX, statistic = mean.boot, R=999)
class(VBLTX.mean.boot)
names(VBLTX.mean.boot)

# print, plot and qqnorm methods
VBLTX.mean.boot
se.muhat["VBLTX"]

# plot bootstrap distribution and qq-plot against normal
plot(VBLTX.mean.boot)


# compute bootstrap confidence intervals from normal approximation
# basic bootstrap method and percentile intervals
boot.ci(VBLTX.mean.boot, conf = 0.95, type = c("norm","perc"))

#
# boostrap SD estimate
#
# function for bootstrapping sample standard deviation
sd.boot = function(x, idx) {
# arguments:
# x 		data to be resampled
# idx		vector of scrambled indices created by boot() function
# value:
# ans		sd value computed using resampled data
     ans = sd(x[idx])
     ans
}

VBLTX.sd.boot = boot(VBLTX, statistic = sd.boot, R=999)
VBLTX.sd.boot
se.sigmahat["VBLTX"]

# plot bootstrap distribution
plot(VBLTX.sd.boot)

# compute confidence intervals
boot.ci(VBLTX.sd.boot, conf=0.95, type=c("norm", "basic", "perc"))

# bootstrap correlation

# function to compute correlation between 1st 2 variables in matrix
rho.boot = function(x.mat, idx) {
# x.mat	n x 2 data matrix to be resampled
# idx		vector of scrambled indices created by boot() function
# value:
# ans		correlation value computed using resampled data

	ans = cor(x.mat[idx,])[1,2]
	ans
}
VBLTX.FMAGX.cor.boot = boot(ret.mat[,c("VBLTX","FMAGX")],
                           statistic=rho.boot, R = 999)
VBLTX.FMAGX.cor.boot
se.rhohat[1]

# plot bootstrap distribution
plot(VBLTX.FMAGX.cor.boot)

# bootstrap confidence intervals
boot.ci(VBLTX.FMAGX.cor.boot, conf=0.95, type=c("norm", "perc"))

#
# 8. Bootstrap VaR
#

# 5% Value-at-Risk
ValueAtRisk.boot = function(x, idx, p=0.05, w=100000) {
# x.mat	data to be resampled
# idx		vector of scrambled indices created by boot() function
# p		probability value for VaR calculation
# w		value of initial investment
# value:
# ans		Value-at-Risk computed using resampled data

	q = mean(x[idx]) + sd(x[idx])*qnorm(p)
	VaR = (exp(q) - 1)*w
	VaR
}

VBLTX.VaR.boot = boot(VBLTX, statistic = ValueAtRisk.boot, R=999)
VBLTX.VaR.boot
boot.ci(VBLTX.VaR.boot, conf=0.95, type=c("norm", "perc"))
plot(VBLTX.VaR.boot)