#
# 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) 

# P1 (Fit beta models to SNP 500 stocks): 
#
dat = read.csv("../../BookCode/Data/Stock_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


# P 7
#
AlphaBeta = read.csv('../../BookCode/Data/AlphaBeta.csv', header=TRUE)
colnames(AlphaBeta) = c('stock', 'alpha', 'beta')
head(AlphaBeta)

#AlphaBeta = AlphaBeta[1:3, ] # for debugging
#AlphaBeta = AlphaBeta[1:20, ] # for Problem 8 

print( c(min(AlphaBeta$beta), max(AlphaBeta$beta)) )

library(linprog)

# Note: B1 is the UPPER bound on the weights in the portfolio and
#      -B2 is the LOWER bound on the weights in the portfolio that is:
#
# -B2 <= wi <= +B1
#
B1 = 0.25
B2 = 0.25

M = dim(AlphaBeta)[1]

AmatLP1 = cbind(diag(1, nrow=M), matrix(0, nrow=M, ncol=M))
AmatLP2 = cbind(matrix(0, nrow=M, ncol=M), diag(1, nrow=M))
AmatLP3 = c(AlphaBeta$beta, -AlphaBeta$beta)
AmatLP = rbind(AmatLP1, AmatLP2, AmatLP3)
bvecLP = c(rep(B1, M), rep(B2, M), 0)
cLP = c(AlphaBeta$alpha, -AlphaBeta$alpha)
const.dir = c(rep('<=', 2*M),'=')

resultLP = solveLP(cvec=cLP, bvec=bvecLP, Amat=AmatLP, lpSolve=TRUE, const.dir=const.dir, maximum=TRUE)

sol = resultLP$solution
x1 = sol[1:M]
x2 = sol[(M+1):(2*M)]
w= x1 - x2 # these are the portfolio weights
print(sprintf('alpha_F (max)= %f; beta_P= %f', resultLP$opt, sum(w * AlphaBeta$beta)))