# mu sigma optimale Portfolio 2D
# Autor: Manfred Jaeger-Ambrozewicz
# Version 20260118

rm(list=ls())

###### Die Daten ########

ones = c(1,1)
mu = c(0.05, 0.01)
sigma = matrix( c(1/4 , 1/10 ,
                  1/10, 4/25), 
                nrow=2, ncol=2)

rf = 0.005
mustern = 0.04

round(sigma,5)
# Standardabweichungen in Prozent
round((sqrt(diag(sigma))),4)*100

###### Sigma invertieren  ########

sigmainv = solve(sigma)
sigmainv
round(sigmainv,2)


###### GMV Portfolio #####

aux1 = 1/as.numeric(t(ones)%*%sigmainv%*%ones)
aux1
sigmainv%*%ones

w1 = aux1*sigmainv%*%ones 
w1

rw1 = t(w1)%*%mu
rw1

sdw1 = sqrt(t(w1)%*%sigma%*%w1)
sdw1

##### Tangentialportfolio ######

aux2 = as.numeric(1/(t(ones)%*%sigmainv%*%(mu-rf*ones)))
aux2

sigmainv%*%(mu-rf*ones)

w2 = aux2*sigmainv%*%(mu-rf*ones)
w2
sum(w2)

rw2 = t(w2)%*%mu
rw2

t(w2)%*%sigma%*%w2

sdw2 = sqrt(t(w2)%*%sigma%*%w2)
sdw2

effw2 = as.numeric(ifelse(rw2 - rf > 0, 1, -1))  
srw2 = (rw2-rf)/sdw2

###### Grenzportfolio mit Zielrendite #####
# ohne die risikolose Anlageform

M = matrix(c(t(mu),t(ones)),ncol=2)
M
B = t(M)%*%sigmainv%*%M
Binv = solve(B)
Binv
M%*%Binv%*%c(mustern,1)

w3 = sigmainv%*%M%*%Binv%*%c(mustern,1)
w3
sum(w3)

rw3 = t(w3)%*%mu 
rw3
sdw3 = sqrt(t(w3)%*%sigma%*%w3)
sdw3

###### Grenzportfolio mit Zielrendite ######
# einschliesslich risikolose Anlageform

aux3 = as.numeric( (mustern-rf)/( t(mu-rf*ones)%*%sigmainv%*%(mu-rf*ones) ) )
aux3

w4 = aux3*sigmainv%*%(mu-rf*ones)
w4


w4s = c(w4,1-sum(w4))
w4s
mus = c(mu, rf)

sum(w4s*mus) 
sum(w4)

sum(w4) 
alpha = (mustern-rf)/(rw2-rf)
alpha

w4alternativ = as.numeric(alpha)*w2
w4alternativ
w4

rw4 = alpha*rw2 + (1-alpha)*rf
rw4

sdw4 = sqrt(alpha^2*sdw2^2)
alpha*sdw2
sdw4

sqrt(t(w4)%*%sigma%*%w4)
t(w4)%*%sigma%*%w4

w4/as.numeric(alpha)
w2
as.numeric(alpha)*w2
w4

##### Abbildungen #####

A = as.numeric(t(mu)%*%sigmainv%*%mu)
B = as.numeric(t(mu)%*%sigmainv%*%ones)
C = as.numeric(t(ones)%*%sigmainv%*%ones)

hyperbel = function(x){
  out = sqrt((A-2*B*x+C*x^2)/(A*C-B^2))
}

weite = max(abs(c(rw1)-mu))*2
schritte = 2*weite/100
Rendite1 = seq(from = c(rw1) - weite, to = rw1, by = schritte )
Standardabweichungen1 = hyperbel(Rendite1)
Rendite2 = seq(from = c(rw1), to = rw1 + weite, by = schritte )
Standardabweichungen2 = hyperbel(Rendite2)
weite2 = 1.3*max(Standardabweichungen1,Standardabweichungen2)

print(hyperbel(0.045))

par(mai = c(1.2,1.2,0.5,1.2))
plot(Standardabweichungen1,Rendite1,
     xlim = c(0,weite2),
     ylim = c(rw1-weite,rw1+weite),
     bty = "L",
     type = "l",
     xlab = "Standardabweichung",
     ylab = "erw. Rendite",
     cex.lab = 1.3,
     col = "red",
     lwd = 3)

points(Standardabweichungen2,Rendite2,
       type = "l",
       col  = "green",
       lwd = 3)
points(sdw1,rw1,
       col = "blue",
       cex = 2,
       pch = 8)
points(sqrt(diag(sigma)),mu,
       col = "yellow",
       cex = 2,
       pch = 16)
points(sdw3,rw3,
       col = "black",
       cex = 2,
       pch = 16)
abline(0,0,col="grey")

weite = max(abs(c(rw1)-mu),rf)*2
schritte = 2*weite/100
Rendite1 = seq(from = c(rw1) - weite, to = rw1, by = schritte )
Standardabweichungen1 = hyperbel(Rendite1)
Rendite2 = seq(from = c(rw1), to = rw1 + weite, by = schritte )
Standardabweichungen2 = hyperbel(Rendite2)
weite2 = 1.3*max(Standardabweichungen1,Standardabweichungen2)

par(mai = c(1.2,1.2,0.5,1.2))
plot(Standardabweichungen1,Rendite1,
     xlim = c(0,weite2),
     ylim = c(rw1-weite,rw1+weite),
     bty = "L",
     type = "l",
     xlab = "Standardabweichung",
     ylab = "erw. Rendite",
     cex.lab = 1.3,
     col = "red",
     lwd = 3)

points(Standardabweichungen2,Rendite2,
       type = "l",
       col  = "green",
       lwd = 3)
points(sdw1,rw1,
       col = "black",
       cex = 2,
       pch = 8)
points(sdw2,rw2,
       col = "blue",
       cex = 2,
       pch = 8)
points(sqrt(diag(sigma)),mu,
       col = "yellow",
       cex = 2,
       pch = 16)
points(sdw2,rw2,
       col = "blue",
       cex = 2,
       pch = 8)
points(sdw3,rw3,
       col = "black",
       cex = 2,
       pch = 16)
points(0,rf,
       col = "blue",
       cex = 2,
       pch = 8)
lines(c(0,Standardabweichungen1),
      c(rf,rf+effw2*c(srw2)*Standardabweichungen1),lwd=2,lty=1,
      col = "green")
lines(c(0,Standardabweichungen1),
      c(rf,rf-effw2*c(srw2)*Standardabweichungen1),lwd=2,lty=2,
      col = "red")
points(sdw4,rw4,
       col = "green",
       cex = 2,
       pch = 8)
abline(0,0,col="grey")