#install.packages("PerformanceAnalytics", repos = "http://cran.us.r-project.org")
#install.packages("dplyr", repos = "http://cran.us.r-project.org")
#install.packages("tidyquant", repos = "http://cran.us.r-project.org")
#install.packages("quantmod", repos = "http://cran.us.r-project.org")
#install.packages("tseries", repos = "http://cran.us.r-project.org")
#install.packages("tidyverse", repos = "http://cran.us.r-project.org")



library(PerformanceAnalytics)  #useful package !!!
library(quantmod)              # useful package
library(tidyquant)             # useful

# other libraries

library(ggplot2)
library(dplyr)

library(tseries)
library(tidyverse)
library(plotly)
library(hrbrthemes)
library(xts)
library(knitr)
library(kableExtra)
library(car)
library(mathjaxr)
library(zoo)
rm(list=ls())
symbol_name <<- c("AAPL", "GOOG", "AMZN", "F", "T", "TQQQ")
#symbol_name <- "AAPL"

for (i in 1:length(symbol_name)) {
  prac <<- Ad(getSymbols(symbol_name[i], from = "2020-01-01", to = "2022-12-31",auto.assign=FALSE))
  if (i==1) {
    price <<-prac
  } else{
    price <<- merge(price,prac)
  }
}
rm(prac)    # prac is just temporary variable to remove
colnames(price) <- symbol_name  #puting the names of the shares
returns <<- CalculateReturns(price, method="log")
for(i in 1:dim(returns)[2]){   #imputation of the missing data
  returns[,i][is.na(returns[,i])] <- median(returns[,i],na.rm = TRUE)
}

Some preliminary calculations

Correlation structure of the portfolio assets

If speaking about facing the risks, the effective diversification of the investment is needed. In this way, we can reduce the unsystemic risks. See the correlations in our hypotetical portfolio in the following Figure.

#chart.Correlation(returns[2:ncol(return_p)])
chart.Correlation(returns)

Estimation of the variance-covariance matrix of the portfolio returns

Variance - covariance matrix is a central topic of the portfolio theory. It containes the return variances on the main diagonal, while besides there are the covariances. If speaking about the correlation matrix, then the diagonal terms are 1’s and out of the diagonal, there are the corresponding correlations.

var_covar_p <- cov(returns)
var_covar_p
             AAPL         GOOG         AMZN            F            T         TQQQ
AAPL 0.0005400903 0.0003680854 0.0003809038 0.0003066046 0.0001699712 0.0011602777
GOOG 0.0003680854 0.0004691521 0.0003656386 0.0002892226 0.0001543824 0.0010460374
AMZN 0.0003809038 0.0003656386 0.0006058288 0.0002275790 0.0001062797 0.0010926341
F    0.0003066046 0.0002892226 0.0002275790 0.0009622201 0.0002537433 0.0008846587
T    0.0001699712 0.0001543824 0.0001062797 0.0002537433 0.0003172865 0.0004525939
TQQQ 0.0011602777 0.0010460374 0.0010926341 0.0008846587 0.0004525939 0.0031640329

Markowitz portfolio

mean_p <- apply(returns, 2, mean)
mean_p
         AAPL          GOOG          AMZN             F             T          TQQQ 
 0.0007532886  0.0003463638 -0.0001606586  0.0003487300 -0.0002794560 -0.0003477078 
# we will draw 100 points to show Markowitz set
NAH <- matrix(runif(10000*(dim(returns)[2]-1)),ncol=(dim(returns)[2]-1))
NAH_sorted <- t(apply(NAH, 1, sort))
MATICA_ALFA <- matrix(rep(0,10000*dim(returns)[2]),ncol=dim(returns)[2])
for (i in 1:10000) {
  MATICA_ALFA[i,1] = NAH_sorted[i,1]
  MATICA_ALFA[i,2] = NAH_sorted[i,2] - NAH_sorted[i,1]
  MATICA_ALFA[i,3] = NAH_sorted[i,3] - NAH_sorted[i,2]
  MATICA_ALFA[i,4] = NAH_sorted[i,4] - NAH_sorted[i,3]
  MATICA_ALFA[i,5] = NAH_sorted[i,5] - NAH_sorted[i,4]
  MATICA_ALFA[i,6] = 1 - NAH_sorted[i,5]
}

MeanReturns <- rep(0,times=10000)
StdReturns <- rep(0,times=10000)
for (i in 1:10000){
  MeanReturns[i] <- MATICA_ALFA[i,]%*%mean_p
  StdReturns[i] <- sqrt(t(MATICA_ALFA[i,])%*%var_covar_p%*%t(t(MATICA_ALFA[i,])))
}

DataForPicture <- data.frame(cbind(StdReturns,MeanReturns))


ggplot(DataForPicture, aes(x=StdReturns, y=MeanReturns)) + 
    geom_point(colour = "red", size = 1,alpha=1/10)

rather simple optimization

returns <<- returns
Optimalne_Portfolio <- portfolio.optim(returns,pm=0.0001,shorts=TRUE)
Optimalne_Portfolio$pw
[1]  0.11799636  0.27237890  0.42328060  0.02451135  0.47616119 -0.31432839
Optimalne_Portfolio$pm
[1] 1e-04
Optimalne_Portfolio$ps
[1] 0.009759287

The constraint for mean value is is 10^{-5} with the standard deviation of 0.0097593 and the optimal composition0.1179964, 0.2723789, 0.4232806, 0.0245113, 0.4761612, -0.3143284.

#returns_matrix <- as.matrix(returns)

#assets <- colnames(returns_matrix)

# Create the portfolio specification
#portfolio_spec <- portfolio.spec(assets)

# Add constraints to prevent short selling
#portfolio_spec <- add.constraint(portfolio_spec, type = "long_only",min_sum=0.99,max_sum=1.01)
#portfolio.spec <- add.constraint(portfolio_spec, type="full_investment")
#portfolio_spec <- add.constraint(portfolio_spec, type="weight_sum", min_sum=1, max_sum=1)
#add.constraint(portfolio_spec, type="weight_sum", min_sum=0.99, max_sum=1.01, indexnum=1)

# Set the objective function to minimize variance
#portfolio_spec <- add.objective(portfolio_spec, type = "risk", name = "var")

# definition of the optimization problem
#portfolio <-optimize.portfolio(returns_matrix, portfolio_spec,optimize_method="random",trace=TRUE)
#portfolio <- optimize.portfolio(returns, portfolio_spec)

# view the portfolio allocation
#print(portfolio)
#chart.EfficientFrontier(portfolio,match.col="StdDev")
LS0tCnRpdGxlOiAiUG9ydGZvbGlvIDMiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCmBgYHtyIGxpYnJhcmllc30KI2luc3RhbGwucGFja2FnZXMoIlBlcmZvcm1hbmNlQW5hbHl0aWNzIiwgcmVwb3MgPSAiaHR0cDovL2NyYW4udXMuci1wcm9qZWN0Lm9yZyIpCiNpbnN0YWxsLnBhY2thZ2VzKCJkcGx5ciIsIHJlcG9zID0gImh0dHA6Ly9jcmFuLnVzLnItcHJvamVjdC5vcmciKQojaW5zdGFsbC5wYWNrYWdlcygidGlkeXF1YW50IiwgcmVwb3MgPSAiaHR0cDovL2NyYW4udXMuci1wcm9qZWN0Lm9yZyIpCiNpbnN0YWxsLnBhY2thZ2VzKCJxdWFudG1vZCIsIHJlcG9zID0gImh0dHA6Ly9jcmFuLnVzLnItcHJvamVjdC5vcmciKQojaW5zdGFsbC5wYWNrYWdlcygidHNlcmllcyIsIHJlcG9zID0gImh0dHA6Ly9jcmFuLnVzLnItcHJvamVjdC5vcmciKQojaW5zdGFsbC5wYWNrYWdlcygidGlkeXZlcnNlIiwgcmVwb3MgPSAiaHR0cDovL2NyYW4udXMuci1wcm9qZWN0Lm9yZyIpCgoKCmxpYnJhcnkoUGVyZm9ybWFuY2VBbmFseXRpY3MpICAjdXNlZnVsIHBhY2thZ2UgISEhCmxpYnJhcnkocXVhbnRtb2QpICAgICAgICAgICAgICAjIHVzZWZ1bCBwYWNrYWdlCmxpYnJhcnkodGlkeXF1YW50KSAgICAgICAgICAgICAjIHVzZWZ1bAoKIyBvdGhlciBsaWJyYXJpZXMKCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShkcGx5cikKCmxpYnJhcnkodHNlcmllcykKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkocGxvdGx5KQpsaWJyYXJ5KGhyYnJ0aGVtZXMpCmxpYnJhcnkoeHRzKQpsaWJyYXJ5KGtuaXRyKQpsaWJyYXJ5KGthYmxlRXh0cmEpCmxpYnJhcnkoY2FyKQpsaWJyYXJ5KG1hdGhqYXhyKQpsaWJyYXJ5KHpvbykKcm0obGlzdD1scygpKQpgYGAKCgoKCmBgYHtyIGRvd25sb2FkX2RhdGEsIG1lc3NhZ2U9RkFMU0V9CnN5bWJvbF9uYW1lIDw8LSBjKCJBQVBMIiwgIkdPT0ciLCAiQU1aTiIsICJGIiwgIlQiLCAiVFFRUSIpCiNzeW1ib2xfbmFtZSA8LSAiQUFQTCIKCmZvciAoaSBpbiAxOmxlbmd0aChzeW1ib2xfbmFtZSkpIHsKICBwcmFjIDw8LSBBZChnZXRTeW1ib2xzKHN5bWJvbF9uYW1lW2ldLCBmcm9tID0gIjIwMjAtMDEtMDEiLCB0byA9ICIyMDIyLTEyLTMxIixhdXRvLmFzc2lnbj1GQUxTRSkpCiAgaWYgKGk9PTEpIHsKICAgIHByaWNlIDw8LXByYWMKICB9IGVsc2V7CiAgICBwcmljZSA8PC0gbWVyZ2UocHJpY2UscHJhYykKICB9Cn0Kcm0ocHJhYykgICAgIyBwcmFjIGlzIGp1c3QgdGVtcG9yYXJ5IHZhcmlhYmxlIHRvIHJlbW92ZQpjb2xuYW1lcyhwcmljZSkgPC0gc3ltYm9sX25hbWUgICNwdXRpbmcgdGhlIG5hbWVzIG9mIHRoZSBzaGFyZXMKYGBgCgoKYGBge3IgQ29tcHV0YXRpb25fb2ZfcmV0dXJuc30KcmV0dXJucyA8PC0gQ2FsY3VsYXRlUmV0dXJucyhwcmljZSwgbWV0aG9kPSJsb2ciKQpmb3IoaSBpbiAxOmRpbShyZXR1cm5zKVsyXSl7ICAgI2ltcHV0YXRpb24gb2YgdGhlIG1pc3NpbmcgZGF0YQogIHJldHVybnNbLGldW2lzLm5hKHJldHVybnNbLGldKV0gPC0gbWVkaWFuKHJldHVybnNbLGldLG5hLnJtID0gVFJVRSkKfQpgYGAKCiMjIFNvbWUgcHJlbGltaW5hcnkgY2FsY3VsYXRpb25zCgoKCgoKIyMgQ29ycmVsYXRpb24gc3RydWN0dXJlIG9mIHRoZSBwb3J0Zm9saW8gYXNzZXRzCgpJZiBzcGVha2luZyBhYm91dCBmYWNpbmcgdGhlIHJpc2tzLCB0aGUgZWZmZWN0aXZlIGRpdmVyc2lmaWNhdGlvbiBvZiB0aGUgaW52ZXN0bWVudCBpcyBuZWVkZWQuIEluIHRoaXMgd2F5LCB3ZSBjYW4gcmVkdWNlIHRoZSAqdW5zeXN0ZW1pYyByaXNrcyouIFNlZSB0aGUgY29ycmVsYXRpb25zIGluIG91ciBoeXBvdGV0aWNhbCBwb3J0Zm9saW8gaW4gdGhlIGZvbGxvd2luZyBGaWd1cmUuCgoKYGBge3IgZmlnNCwgZXJyb3I9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CiNjaGFydC5Db3JyZWxhdGlvbihyZXR1cm5zWzI6bmNvbChyZXR1cm5fcCldKQpjaGFydC5Db3JyZWxhdGlvbihyZXR1cm5zKQpgYGAKIyMgRXN0aW1hdGlvbiBvZiB0aGUgdmFyaWFuY2UtY292YXJpYW5jZSBtYXRyaXggb2YgdGhlIHBvcnRmb2xpbyByZXR1cm5zCgpWYXJpYW5jZSAtIGNvdmFyaWFuY2UgbWF0cml4IGlzIGEgY2VudHJhbCB0b3BpYyBvZiB0aGUgcG9ydGZvbGlvIHRoZW9yeS4gSXQgY29udGFpbmVzIHRoZSByZXR1cm4gdmFyaWFuY2VzIG9uIHRoZSBtYWluIGRpYWdvbmFsLCB3aGlsZSBiZXNpZGVzIHRoZXJlIGFyZSB0aGUgY292YXJpYW5jZXMuIElmIHNwZWFraW5nIGFib3V0IHRoZSBjb3JyZWxhdGlvbiBtYXRyaXgsIHRoZW4gdGhlIGRpYWdvbmFsIHRlcm1zIGFyZSAxJ3MgYW5kIG91dCBvZiB0aGUgZGlhZ29uYWwsIHRoZXJlIGFyZSB0aGUgY29ycmVzcG9uZGluZyBjb3JyZWxhdGlvbnMuCgpgYGB7ciB2YXJpYW5jZS1jb3ZhcmlhbmNlX21hdHJpeH0KdmFyX2NvdmFyX3AgPC0gY292KHJldHVybnMpCnZhcl9jb3Zhcl9wCmBgYAoKIyMgTWFya293aXR6IHBvcnRmb2xpbwoKYGBge3J9Cm1lYW5fcCA8LSBhcHBseShyZXR1cm5zLCAyLCBtZWFuKQptZWFuX3AKYGBgCgpgYGB7cn0KIyB3ZSB3aWxsIGRyYXcgMTAwIHBvaW50cyB0byBzaG93IE1hcmtvd2l0eiBzZXQKTkFIIDwtIG1hdHJpeChydW5pZigxMDAwMCooZGltKHJldHVybnMpWzJdLTEpKSxuY29sPShkaW0ocmV0dXJucylbMl0tMSkpCk5BSF9zb3J0ZWQgPC0gdChhcHBseShOQUgsIDEsIHNvcnQpKQpNQVRJQ0FfQUxGQSA8LSBtYXRyaXgocmVwKDAsMTAwMDAqZGltKHJldHVybnMpWzJdKSxuY29sPWRpbShyZXR1cm5zKVsyXSkKZm9yIChpIGluIDE6MTAwMDApIHsKICBNQVRJQ0FfQUxGQVtpLDFdID0gTkFIX3NvcnRlZFtpLDFdCiAgTUFUSUNBX0FMRkFbaSwyXSA9IE5BSF9zb3J0ZWRbaSwyXSAtIE5BSF9zb3J0ZWRbaSwxXQogIE1BVElDQV9BTEZBW2ksM10gPSBOQUhfc29ydGVkW2ksM10gLSBOQUhfc29ydGVkW2ksMl0KICBNQVRJQ0FfQUxGQVtpLDRdID0gTkFIX3NvcnRlZFtpLDRdIC0gTkFIX3NvcnRlZFtpLDNdCiAgTUFUSUNBX0FMRkFbaSw1XSA9IE5BSF9zb3J0ZWRbaSw1XSAtIE5BSF9zb3J0ZWRbaSw0XQogIE1BVElDQV9BTEZBW2ksNl0gPSAxIC0gTkFIX3NvcnRlZFtpLDVdCn0KCk1lYW5SZXR1cm5zIDwtIHJlcCgwLHRpbWVzPTEwMDAwKQpTdGRSZXR1cm5zIDwtIHJlcCgwLHRpbWVzPTEwMDAwKQpmb3IgKGkgaW4gMToxMDAwMCl7CiAgTWVhblJldHVybnNbaV0gPC0gTUFUSUNBX0FMRkFbaSxdJSolbWVhbl9wCiAgU3RkUmV0dXJuc1tpXSA8LSBzcXJ0KHQoTUFUSUNBX0FMRkFbaSxdKSUqJXZhcl9jb3Zhcl9wJSoldCh0KE1BVElDQV9BTEZBW2ksXSkpKQp9CgpEYXRhRm9yUGljdHVyZSA8LSBkYXRhLmZyYW1lKGNiaW5kKFN0ZFJldHVybnMsTWVhblJldHVybnMpKQoKCmdncGxvdChEYXRhRm9yUGljdHVyZSwgYWVzKHg9U3RkUmV0dXJucywgeT1NZWFuUmV0dXJucykpICsgCiAgICBnZW9tX3BvaW50KGNvbG91ciA9ICJyZWQiLCBzaXplID0gMSxhbHBoYT0xLzEwKQoKYGBgCgoKCgoKCgojIyByYXRoZXIgc2ltcGxlIG9wdGltaXphdGlvbgoKYGBge3J9CnJldHVybnMgPDwtIHJldHVybnMKT3B0aW1hbG5lX1BvcnRmb2xpbyA8LSBwb3J0Zm9saW8ub3B0aW0ocmV0dXJucyxwbT0wLjAwMDEsc2hvcnRzPVRSVUUpCk9wdGltYWxuZV9Qb3J0Zm9saW8kcHcKT3B0aW1hbG5lX1BvcnRmb2xpbyRwbQpPcHRpbWFsbmVfUG9ydGZvbGlvJHBzCmBgYApUaGUgY29uc3RyYWludCBmb3IgbWVhbiB2YWx1ZSBpcyBpcyBgciBPcHRpbWFsbmVfUG9ydGZvbGlvJHBtYCB3aXRoIHRoZSBzdGFuZGFyZCBkZXZpYXRpb24gb2YgYHIgT3B0aW1hbG5lX1BvcnRmb2xpbyRwc2AgYW5kIHRoZSBvcHRpbWFsIGNvbXBvc2l0aW9uYHIgT3B0aW1hbG5lX1BvcnRmb2xpbyRwd2AuCgoKCgoKCmBgYHtyIHJvenJvYmVuZX0KI3JldHVybnNfbWF0cml4IDwtIGFzLm1hdHJpeChyZXR1cm5zKQoKI2Fzc2V0cyA8LSBjb2xuYW1lcyhyZXR1cm5zX21hdHJpeCkKCiMgQ3JlYXRlIHRoZSBwb3J0Zm9saW8gc3BlY2lmaWNhdGlvbgojcG9ydGZvbGlvX3NwZWMgPC0gcG9ydGZvbGlvLnNwZWMoYXNzZXRzKQoKIyBBZGQgY29uc3RyYWludHMgdG8gcHJldmVudCBzaG9ydCBzZWxsaW5nCiNwb3J0Zm9saW9fc3BlYyA8LSBhZGQuY29uc3RyYWludChwb3J0Zm9saW9fc3BlYywgdHlwZSA9ICJsb25nX29ubHkiLG1pbl9zdW09MC45OSxtYXhfc3VtPTEuMDEpCiNwb3J0Zm9saW8uc3BlYyA8LSBhZGQuY29uc3RyYWludChwb3J0Zm9saW9fc3BlYywgdHlwZT0iZnVsbF9pbnZlc3RtZW50IikKI3BvcnRmb2xpb19zcGVjIDwtIGFkZC5jb25zdHJhaW50KHBvcnRmb2xpb19zcGVjLCB0eXBlPSJ3ZWlnaHRfc3VtIiwgbWluX3N1bT0xLCBtYXhfc3VtPTEpCiNhZGQuY29uc3RyYWludChwb3J0Zm9saW9fc3BlYywgdHlwZT0id2VpZ2h0X3N1bSIsIG1pbl9zdW09MC45OSwgbWF4X3N1bT0xLjAxLCBpbmRleG51bT0xKQoKIyBTZXQgdGhlIG9iamVjdGl2ZSBmdW5jdGlvbiB0byBtaW5pbWl6ZSB2YXJpYW5jZQojcG9ydGZvbGlvX3NwZWMgPC0gYWRkLm9iamVjdGl2ZShwb3J0Zm9saW9fc3BlYywgdHlwZSA9ICJyaXNrIiwgbmFtZSA9ICJ2YXIiKQoKIyBkZWZpbml0aW9uIG9mIHRoZSBvcHRpbWl6YXRpb24gcHJvYmxlbQojcG9ydGZvbGlvIDwtb3B0aW1pemUucG9ydGZvbGlvKHJldHVybnNfbWF0cml4LCBwb3J0Zm9saW9fc3BlYyxvcHRpbWl6ZV9tZXRob2Q9InJhbmRvbSIsdHJhY2U9VFJVRSkKI3BvcnRmb2xpbyA8LSBvcHRpbWl6ZS5wb3J0Zm9saW8ocmV0dXJucywgcG9ydGZvbGlvX3NwZWMpCgojIHZpZXcgdGhlIHBvcnRmb2xpbyBhbGxvY2F0aW9uCiNwcmludChwb3J0Zm9saW8pCiNjaGFydC5FZmZpY2llbnRGcm9udGllcihwb3J0Zm9saW8sbWF0Y2guY29sPSJTdGREZXYiKQoKCmBgYAoK