#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

save_plots = F

set.seed(0) 

#
# Fit beta models to SNP 500
#
dat = read.csv("../../BookCode/Data/Stock_FX_Bond_2004_to_2006.csv",header=T)
prices = dat[,c(5,7,9,11,13,15,17,24)]
n = dim(prices)[1] 

dat2 = as.matrix( cbind( dat[(2:n),3]/365, 100*(prices[2:n,]/prices[1:(n-1),] - 1) ) )
colnames(dat2)[1] = "treasury" 
risk_free = dat2[,1]
ExRet = dat2[,2:9] - risk_free
market = ExRet[,8]
stockExRet = ExRet[,1:7]

fit_reg = lm( stockExRet ~ market ) 
summary(fit_reg)
res = residuals(fit_reg)
pairs(res)
options(digits=3)
betas = fit_reg$coeff[2,] 

# P2 EPage 439
#
# Eq. 16.19 is
#
# R_{j,t} = mu_{f,t} + beta_j ( R_{M,t} - mu_{f,t} ) + epsilon_{j,t} 
#
mean(market) * betas

# vs.
#
apply( stockExRet, 2, mean )


# P3 EPage 439
#
cor(res)

# P4
#
summary_fit_reg = summary(fit_reg) 
sigmas = sapply( summary_fit_reg, function(x){ x$sigma } ) # Extract the idiosyncratic risk
names(sigmas) = names(betas)
sigmas2 = sigmas^2
sigma2_M = var(market) 
as.matrix(betas) %*% t(as.matrix(betas)) * sigma2_M + diag( sigmas2 )

# vs.
#
cov( stockExRet )


# P 5
#
summary_fit_reg[3]


# P 6
#
ex_ret_M_predicted = 0.04
betas * ex_ret_M_predicted