#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