## read raw data
## Markowitz's MVO model - 8.1.1 Example
## US stocks, bonds and cash: S&P 500 index, 10-year Treasury bond index, 1-day federal fund rate
y <- read.table("chapter8_example1.dat", header = TRUE)

## transform to ts
mvo <- ts(as.matrix(y[,-1]), start = c(1960,1), freq = 1)

## compute returns (r_t - rt-1)/rt-1
mvo.r <- (mvo[2:nrow(mvo),] - mvo[1:nrow(mvo)-1,])/mvo[1:nrow(mvo)-1,]

## arithmetric mean

colMeans(mvo.r)

## geometric mean
geo.mean <- function (x) {
            y<-vector(length=ncol(x))
            for (i in 1:ncol(x)){
            y[i] <- prod(x[,i]+1)^(1/nrow(x))-1
            }
            return(y)
	    }
geo.mean(mvo.r)


exp(mean(log(mvo.r)))

## covariance matrix
cov(mvo.r)

## volatility
sd(mvo.r)

## correlation
cor(mvo.r)

## QP portfolio optimization
library("ROI.obsolete")

foo <- Q_objective(Q = cov(mvo.r), L = rep(0, ncol(mvo.r)))
r<-0.065
full_invest <- L_constraint(rep(1, ncol(mvo.r)), "==", 1)
target_return <- L_constraint(apply(mvo.r, 2, mean), "==", r)
op <- QP(objective = foo, constraints = rbind(full_invest, target_return))
op
sol <- ROI_solve(op, solver = "quadprog")
w <- sol$solution
round(w, 4)

QP.optimization.1 <- function(dat,R) {
    foo <-  Q_objective(Q = cov(dat), L = rep(0, ncol(dat)))
    w<-matrix(nrow=length(R),ncol=(ncol(dat)+1))
    for (i in 1:length(R)){
    full_invest <- L_constraint(rep(1, ncol(dat)), "==", 1)
    target_return <- L_constraint(apply(dat, 2, mean), "==", R[i])
    op <- QP(objective = foo, constraints = rbind(full_invest, target_return))
    sol <- ROI_solve(op, solver = "quadprog")
    w[i,1] <- R[i]
    w[i,2:(ncol(dat)+1)] <- round(sol$solution,4)
    }
    w
    }

r<-seq(from = 0.065, to = 0.105, by = 0.005)
QP.optimization.1(mvo.r,r)

QP.optimization.2 <- function(dat,R) {
    foo <- Q_objective(Q = cov(dat), L = rep(0, ncol(dat)))
    w<-matrix(nrow=length(R),ncol=(ncol(dat)+1))
    for (i in 1:length(R)){
    full_invest <- L_constraint(rep(1, ncol(dat)), "==", 1)
    target_return <- L_constraint(apply(dat, 2, mean), ">=", R[i])
    biggerzero <- L_constraint(diag(1, nrow=ncol(dat)), rep(">=", ncol(dat)), matrix(0,ncol(dat),1))
    op <- QP(objective = foo, constraints = rbind(full_invest, target_return, biggerzero))
    sol <- ROI_solve(op, solver = "quadprog")
    w[i,1] <- R[i]
    w[i,2:(ncol(dat)+1)] <- round(sol$solution,4)
    }
    w
    }

r<-seq(from = 0.065, to = 0.105, by = 0.005)
qp.sol2<-QP.optimization.2(mvo.r,r)
var.mvo.r<-vector(length=9)
for (i in 1:9)
var.mvo.r[i] <- var(mvo.r%*%qp.sol2[i,2:4])

graph1<-cbind(var.mvo.r*100,r*100)
plot(graph1,xlab="SD %", ylab="Exp.Return %")

matplot(qp.sol2[,2:4]*100, type="l", xlab="Exp.Return %", ylab="Percentage invested")

## Exercise 39

QP.optimization.3 <- function(dat,R) {
    foo <- Q_objective(Q = cov(dat), L = rep(0, ncol(dat)))
    w<-matrix(nrow=length(r),ncol=(ncol(dat)+1))
    for (i in 1:length(R)){
    full_invest <- L_constraint(rep(1, ncol(dat)), "==", 1)
    target_return <- L_constraint(geo.mean(dat), ">=", r[i])
    biggerzero <- L_constraint(diag(1, nrow=ncol(dat)), rep(">=", ncol(dat)), matrix(0,ncol(dat),1))
    op <- QP(objective = foo, constraints = rbind(full_invest, target_return, biggerzero))
    sol <- ROI_solve(op, solver = "quadprog")
    w[i,1] <- R[i]
    w[i,2:(ncol(dat)+1)] <- round(sol$solution,4)
    }
    w
    }

r<-seq(from = 0.065, to = 0.105, by = 0.005)
qp.sol3<-QP.optimization.3(mvo.r,r)
var.mvo.r<-vector(length=9)
for (i in 1:9)
var.mvo.r[i] <- var(mvo.r%*%qp.sol3[i,2:4])

graph1<-cbind(var.mvo.r*100,r*100)
plot(graph1,xlab="SD %", ylab="Exp.Return %")

matplot(qp.sol3[,2:4]*100, type="l", xlab="Exp.Return %", ylab="Percentage invested")

## Exercise 40 - introducing a 4th asset: NASDAQ Composite Index

y2 <- read.table("chapter8_data02.dat", header = TRUE)
mvo2 <- ts(as.matrix(y2[,-1]), start = c(1960,1), freq = 1)

## compute returns (r_t - rt-1)/rt-1
mvo.r2.h <- (mvo2[2:nrow(mvo2),] - mvo2[1:nrow(mvo2)-1,])/mvo2[1:nrow(mvo2)-1,]

mvo.r2<-cbind(mvo.r,mvo.r2.h)
qp.sol1.2<-QP.optimization.1(mvo.r2,r)
qp.sol2.2<-QP.optimization.2(mvo.r2,r)
qp.sol3.2<-QP.optimization.3(mvo.r2,r)

matplot(qp.sol1.2[,2:5]*100, type="l", xlab="Exp.Return %", ylab="Percentage invested")
matplot(qp.sol2.2[,2:5]*100, type="l", xlab="Exp.Return %", ylab="Percentage invested")
matplot(qp.sol3.2[,2:5]*100, type="l", xlab="Exp.Return %", ylab="Percentage invested")

## Exercise 42 
## Using historical returns of the stocks in the DJIA, estimate
## their mean ši and covariance matrix. Let R be the median of the šis.
## (i) Solve Markowitz's MVO model to construct a portfolio of stocks from
## the DJIA that has expected return at least R.
## (ii) Generate a random value uniformly in the interval [0.95 mu_i; 1.05 mu_i],
## for each stock i. Resolve Markowitz's MVO model with these mean returns,
## instead of šis as in (i). Compare the results obtained in (i) and (ii).
## (iii) Repeat three more times and average the ¯ve portfolios found in (i),
## (ii) and (iii). Compare this portfolio with the one found in (i)

mvo3.h <- EuStockMarkets #daily
nrow(mvo3.h)
h<-seq(1,nrow(mvo3.h),20)
mvo3<-mvo3.h[h,]


mvo.r3 <- (mvo3[2:nrow(mvo3),] - mvo3[1:nrow(mvo3)-1,])/mvo3[1:nrow(mvo3)-1,]
colMeans(mvo.r3)
r<-0.02

QP.optimization.1(mvo.r3,r)
colnames(mvo.r3,r)

r.u <- vector(length=ncol(mvo.r3))
for (i in 1:ncol(mvo.r3))
r.u[i] <- runif(1, mean(mvo.r3[,i]))
r.u

QP.optimization.4 <- function(dat,R,u) {
    foo <- Q_objective(Q = cov(dat), L = rep(0, ncol(dat)))
    w<-matrix(nrow=length(R),ncol=(ncol(dat)+1))
    for (i in 1:length(R)){
    full_invest <- L_constraint(rep(1, ncol(dat)), "==", 1)
    target_return <- L_constraint(u, "==", R[i])
    op <- QP(objective = foo, constraints = rbind(full_invest, target_return))
    sol <- ROI_solve(op, solver = "quadprog")
    w[i,1] <- R[i]
    w[i,2:(ncol(dat)+1)] <- round(sol$solution,4)
    }
    w
    }

QP.optimization.4(mvo.r3,r,r.u)

r.u2 <- matrix(ncol=ncol(mvo.r3),nrow=5)
for (i in 1:5){
    for (j in 1:4){
    r.u2[i,j] <- runif(1, mean(mvo.r3[,j]))}
    }
r.u2

r.u.mean<-matrix(0,5,5)

for (i in 1:5){
    r.u.mean[i,]<-QP.optimization.4(mvo.r3,r,r.u2[i,])
    }

r.u.mean

colMeans(r.u.mean[,2:5])
sum(colMeans(r.u.mean[,2:5]))

## Example: Black-Litterman
## assumption 1: Money Market rate will be 2 % next year
## assumption 2: S&P 500 will outperform 10-year Treasury Bonds by 5 %

t<-0.1
Q<-cov(mvo.r)
P<-matrix(c(0,1,0,-1,1,0),2,3)
O<-matrix(c(0.00001,0,0,0.001),2,2)
q<-matrix(c(0.02,0.05),2,1)
p <- colMeans(mvo.r)

mu <- solve(solve(t*Q)+t(P)%*%solve(O)%*%P)%*%(solve(t*Q)%*%p+t(P)%*%solve(O)%*%t(q))
mu<-t(mu)
r<-seq(from = 0.04, to = 0.115, by = 0.005)

QP.optimization.4(mvo.r,r,mu)
colnames(mvo.r)

qp.sol4<-QP.optimization.4(mvo.r,r,mu)
var.mvo.r4<-vector(length=length(r))
for (i in 1:length(r))
var.mvo.r4[i] <- var(mvo.r%*%qp.sol4[i,2:4])

graph4<-cbind(var.mvo.r4*100,r*100)
plot(graph4,xlab="SD %", ylab="Exp.Return %")

matplot(qp.sol4[,2:4]*100, type="l", xlab="Exp.Return %", ylab="Percentage invested")

## Konno and Yamazaki: Mean-Absolute-Deviation (MAD)
## Large scale portfolio optimization

r<-seq(from = 0.065, to = 0.105, by = 0.005)
??? wie soll das gehen???