1 Data Preparation

rm(list = ls())
library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.1.1     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(tidyquant)

## Loading required package: lubridate

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

## Loading required package: PerformanceAnalytics

## Loading required package: xts

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
## Attaching package: 'xts'

## The following objects are masked from 'package:dplyr':
## 
##     first, last

## 
## Attaching package: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':
## 
##     legend

## Loading required package: quantmod

## Loading required package: TTR

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## == Need to Learn tidyquant? ====================================================
## Business Science offers a 1-hour course - Learning Lab #9: Performance Analysis & Portfolio Optimization with tidyquant!
## </> Learn more at: https://university.business-science.io/p/learning-labs-pro </>

symbols <- c("FB", "AMZN", "NFLX", "GOOG", "^IXIC")

tq_get(x = symbols, from = "2017-01-01", to = "2020-12-31") -> raw_data

# Using pivot_wider function to transform data with Symbol columns and value is price by date

raw_data %>% 
  select(date, symbol, adjusted) %>% 
  pivot_wider(names_from = symbol, values_from = adjusted, values_fill = 0) %>% 
  slice(-which.min(date)) -> data_wide_price

# ymd function to transform 'date' into Date class

date_ymd <- data_wide_price$date
market_data <- data_wide_price$`^IXIC`
stocks_data <- data_wide_price %>% select(symbols[1:4])

# group by 'Symbol' and perform log return transformation

raw_data %>%
  filter(symbol %in% symbols[1:4]) %>% 
  group_by(symbol) %>%
  tq_transmute(select = adjusted, 
               mutate_fun = periodReturn, 
               period = "daily", 
               type = "log",
               col_rename = "daily_return") %>% 
  pivot_wider(names_from = symbol, values_from = daily_return) %>% 
  slice(-which.min(date)) -> data_wide_return

2 Bayesian Approach for Portfolio Optimization Problem

# Remove date column: 

data_wide_return %>% select(-date) -> return_wide_train

# Calculate stock mean: 
mean_ret <- colMeans(return_wide_train) 

# Covariance matrix: 

cov_mat <- cov(return_wide_train)

# Number of observations: 
T_rows <- nrow(return_wide_train)

# Number of stock symbols: 
N <- ncol(return_wide_train)

# Risk-free rate
Rf <- 1.54/100

# Defining the Sharpe ratio function

sharpe_func <- function(w1,w2,w3,w4){
  weight <- c(w1,w2,w3,w4)
  
  # Sharpe ratio numerator (Difference between return of portfolio and return of the market)
  numerator <- numeric()
  
  for (i in 1:N) {
    numerator[i] <- weight[i]*mean_ret[i]
  }
  
  # Sharpe ratio denominator
  
  denominator<- numeric()
  for (i in 1:N) {
    denominator1 <- numeric()
    for (j in 1:N) {
      denominator1[j] <- weight[i]*weight[j]*cov_mat[i,j]
    }
    denominator[i] <-sum(denominator1)
  }

sharpe_func<- (sum(numerator) - Rf)/sum(denominator)
sharpe_func <- sharpe_func - 1e9 * (round(sum(weight),10)-1)^2

 result <- list(Score = sharpe_func, Pred = 0)
 return(result)
 
}

2.1 Bayesian Optimization function

library(rBayesianOptimization)
# Defining search boundary

search_bound <- list(w1 = c(0,1), w2 = c(0,1),
                     w3 = c(0,1), w4 = c(0,1))

# Defining the initial weight 

set.seed(29)
search_grid <- data.frame(w1 = runif(20,0,1), 
                          w2 = runif(20,0,1),
                          w3 = runif(20,0,1),
                          w4 = runif(20,0,1))

bayes_finance_ucb <- BayesianOptimization(FUN = sharpe_func, bounds = search_bound, 
                     init_grid_dt = search_grid, init_points = 0, 
                     n_iter = 10, acq = "ucb")

## elapsed = 0.01   Round = 1   w1 = 0.09968104 w2 = 0.6658932  w3 = 0.6111826  w4 = 0.4156561  Value = -627918305.7817 
## elapsed = 0.00   Round = 2   w1 = 0.2409021  w2 = 0.618313   w3 = 0.002232353    w4 = 0.1553593  Value = -282512.3176 
## elapsed = 0.00   Round = 3   w1 = 0.1032138  w2 = 0.9728763  w3 = 0.9860362  w4 = 0.5596409  Value = -2630128926.6352 
## elapsed = 0.00   Round = 4   w1 = 0.3255839  w2 = 0.9823095  w3 = 0.6829247  w4 = 0.3582216  Value = -1819907846.1318 
## elapsed = 0.00   Round = 5   w1 = 0.5849888  w2 = 0.6216898  w3 = 0.07372197 w4 = 0.5145345  Value = -631921867.3780 
## elapsed = 0.00   Round = 6   w1 = 0.09308678 w2 = 0.3513491  w3 = 0.8676589  w4 = 0.7653593  Value = -1160907474.5325 
## elapsed = 0.00   Round = 7   w1 = 0.8281676  w2 = 0.4874011  w3 = 0.9599468  w4 = 0.5132438  Value = -3199659798.0306 
## elapsed = 0.00   Round = 8   w1 = 0.8663051  w2 = 0.8992033  w3 = 0.7208769  w4 = 0.6422553  Value = -4531110463.8855 
## elapsed = 0.00   Round = 9   w1 = 0.1199924  w2 = 0.6815601  w3 = 0.5997722  w4 = 0.4208255  Value = -675931096.3426 
## elapsed = 0.00   Round = 10  w1 = 0.2331467  w2 = 0.1606649  w3 = 0.6492987  w4 = 0.2934759  Value = -113290317.9838 
## elapsed = 0.00   Round = 11  w1 = 0.9825898  w2 = 0.366115   w3 = 0.4626285  w4 = 0.1177303  Value = -863159427.6984 
## elapsed = 0.00   Round = 12  w1 = 0.3930963  w2 = 0.1318494  w3 = 0.7139501  w4 = 0.892861   Value = -1280873498.9328 
## elapsed = 0.00   Round = 13  w1 = 0.3010381  w2 = 0.07464958 w3 = 0.2014642  w4 = 0.2470552  Value = -30903213.6227 
## elapsed = 0.00   Round = 14  w1 = 0.6314668  w2 = 0.7221163  w3 = 0.621876   w4 = 0.07783208 Value = -1109422245.1608 
## elapsed = 0.00   Round = 15  w1 = 0.1761709  w2 = 0.4922448  w3 = 0.5012273  w4 = 0.02612822 Value = -38326395.3646 
## elapsed = 0.00   Round = 16  w1 = 0.8274101  w2 = 0.09952085 w3 = 0.4914928  w4 = 0.2887795  Value = -500136442.7797 
## elapsed = 0.00   Round = 17  w1 = 0.6625275  w2 = 0.9217126  w3 = 0.5167407  w4 = 0.8971432  Value = -3992499687.3719 
## elapsed = 0.00   Round = 18  w1 = 0.3621201  w2 = 0.315932   w3 = 0.8234472  w4 = 0.05974056 Value = -314990198.3589 
## elapsed = 0.00   Round = 19  w1 = 0.8667453  w2 = 0.3276701  w3 = 0.5613287  w4 = 0.007461678    Value = -582483115.7811 
## elapsed = 0.00   Round = 20  w1 = 0.3626263  w2 = 0.06703095 w3 = 0.1473852  w4 = 0.2175045  Value = -42211025.7663 
## elapsed = 0.00   Round = 21  w1 = 2.220446e-16   w2 = 2.220446e-16   w3 = 2.220446e-16   w4 = 1.0000 Value = -47.62243 
## elapsed = 0.00   Round = 22  w1 = 2.220446e-16   w2 = 1.0000 w3 = 2.220446e-16   w4 = 2.220446e-16   Value = -37.67061 
## elapsed = 0.00   Round = 23  w1 = 2.220446e-16   w2 = 2.220446e-16   w3 = 1.0000 w4 = 2.220446e-16   Value = -22.62869 
## elapsed = 0.00   Round = 24  w1 = 1.0000 w2 = 2.220446e-16   w3 = 2.220446e-16   w4 = 2.220446e-16   Value = -31.04747 
## elapsed = 0.00   Round = 25  w1 = 2.220446e-16   w2 = 0.4484008  w3 = 2.220446e-16   w4 = 0.651307   Value = -9941690.1670 
## elapsed = 0.00   Round = 26  w1 = 0.4100879  w2 = 2.220446e-16   w3 = 2.220446e-16   w4 = 0.6557107  Value = -4329489.8845 
## elapsed = 0.00   Round = 27  w1 = 2.220446e-16   w2 = 2.220446e-16   w3 = 0.3679016  w4 = 0.6607144  Value = -818918.7523 
## elapsed = 0.00   Round = 28  w1 = 0.6919332  w2 = 0.3270108  w3 = 2.220446e-16   w4 = 2.220446e-16   Value = -358913.3609 
## elapsed = 0.00   Round = 29  w1 = 0.2199992  w2 = 2.220446e-16   w3 = 0.7676581  w4 = 2.220446e-16   Value = -152370.8239 
## elapsed = 0.00   Round = 30  w1 = 2.964779e-09   w2 = 0.698177   w3 = 0.2267262  w4 = 0.06501359 Value = -101711.8503 
## 
##  Best Parameters Found: 
## Round = 23   w1 = 2.220446e-16   w2 = 2.220446e-16   w3 = 1.0000 w4 = 2.220446e-16   Value = -22.62869

2.2 Bayesian Framework for Mean-Covariance model

# Setting hyperparameters for Bayesian portfolio selection

nu <- N+2
tau <- 250 
stan_parameter<- list(
  T_rows = T_rows,
  N = N,
  nu = nu,
  tau = tau,
  eta = mean_ret,
  Ret = as.matrix(return_wide_train),
  omega = cov_mat*(nu - N - 1)
)

#-------------------------------
# Defining model blocks in Stan 
#-------------------------------

# Data block reads external information 

bayes.stan = 
  "
data {
  // 5 is arbitrary
  int<lower = 5> T_rows;
  int<lower = 2> N;
  real<lower = N - 1> nu;
  real<lower = 0> tau;
  vector[N] eta;
  matrix[T_rows, N] Ret;
  cov_matrix[N] omega;
}

// Parameters block defines the sampling space

parameters {
  vector[N] mu;
  cov_matrix[N] sigma;
}

// Transformed parameters block allows for parameter processing before the posterior distribution is computed 

transformed parameters{
  cov_matrix[N] sigma_scaled;
  sigma_scaled = (1 / tau) * sigma;
}

// Model block defines the posterior distribution

model {
  
  target += inv_wishart_lpdf(sigma | nu, omega);
  target += multi_normal_lpdf(mu | eta, sigma_scaled);
  for(t in 1:T_rows){
    target += multi_normal_lpdf(Ret[t]| mu, sigma);
  }
}

"

# Fitting model with stan using No-U-Turn sampling method (variant of Hamilton MCMC method)

library(rstan)

## Loading required package: StanHeaders

## rstan (Version 2.21.3, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)

## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file

## 
## Attaching package: 'rstan'

## The following object is masked from 'package:tidyr':
## 
##     extract

options(mc.cores = parallel::detectCores()) # Allowing for parallel computing 

rstan_options(auto_write = TRUE)

# Setting the number of markov chains:

chains <- N+2 

# Settting number of iterations each chain:

iter <- 2000

# Using No-U-Turn sampling method:

algorithm<- "NUTS"

# Fitting model
system.time(fitting <- stan(model_code = bayes.stan,
    data = stan_parameter,
    algorithm = algorithm, 
    chains = chains,
    warmup = iter / 2,
    iter = iter,
    seed = 42, 
    verbose = TRUE
  )
)

## 
## TRANSLATING MODEL '2af95bb1eb403f4e3d61982a99281fba' FROM Stan CODE TO C++ CODE NOW.
## successful in parsing the Stan model '2af95bb1eb403f4e3d61982a99281fba'.
## OS: x86_64, mingw32; rstan: 2.21.3; Rcpp: 1.0.8; inline: 0.3.19 
##  >> setting environment variables: 
## LOCAL_LIBS =  "D:/R-4.0.2/R-4.1.2/library/rstan/lib/x64/libStanServices.a" -L"D:/R-4.0.2/R-4.1.2/library/StanHeaders/libs/x64" -lStanHeaders -L"D:/R-4.0.2/R-4.1.2/library/RcppParallel/lib/x64" -ltbb
## PKG_CPPFLAGS =   -I"D:/R-4.0.2/R-4.1.2/library/Rcpp/include/"  -I"D:/R-4.0.2/R-4.1.2/library/RcppEigen/include/"  -I"D:/R-4.0.2/R-4.1.2/library/RcppEigen/include/unsupported"  -I"D:/R-4.0.2/R-4.1.2/library/BH/include" -I"D:/R-4.0.2/R-4.1.2/library/StanHeaders/include/src/"  -I"D:/R-4.0.2/R-4.1.2/library/StanHeaders/include/"  -I"D:/R-4.0.2/R-4.1.2/library/RcppParallel/include/"  -I"D:/R-4.0.2/R-4.1.2/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include "D:/R-4.0.2/R-4.1.2/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp"  -std=c++1y
##  >> Program source :
## 
##    1 : 
##    2 : // includes from the plugin
##    3 : // [[Rcpp::plugins(cpp14)]]
##    4 : 
##    5 : 
##    6 : // user includes
##    7 : #include <Rcpp.h>
##    8 : #include <rstan/io/rlist_ref_var_context.hpp>
##    9 : #include <rstan/io/r_ostream.hpp>
##   10 : #include <rstan/stan_args.hpp>
##   11 : #include <boost/integer/integer_log2.hpp>
##   12 : // Code generated by Stan version 2.21.0
##   13 : 
##   14 : #include <stan/model/model_header.hpp>
##   15 : 
##   16 : namespace model99440854acd_2af95bb1eb403f4e3d61982a99281fba_namespace {
##   17 : 
##   18 : using std::istream;
##   19 : using std::string;
##   20 : using std::stringstream;
##   21 : using std::vector;
##   22 : using stan::io::dump;
##   23 : using stan::math::lgamma;
##   24 : using stan::model::prob_grad;
##   25 : using namespace stan::math;
##   26 : 
##   27 : static int current_statement_begin__;
##   28 : 
##   29 : stan::io::program_reader prog_reader__() {
##   30 :     stan::io::program_reader reader;
##   31 :     reader.add_event(0, 0, "start", "model99440854acd_2af95bb1eb403f4e3d61982a99281fba");
##   32 :     reader.add_event(39, 37, "end", "model99440854acd_2af95bb1eb403f4e3d61982a99281fba");
##   33 :     return reader;
##   34 : }
##   35 : 
##   36 : class model99440854acd_2af95bb1eb403f4e3d61982a99281fba
##   37 :   : public stan::model::model_base_crtp<model99440854acd_2af95bb1eb403f4e3d61982a99281fba> {
##   38 : private:
##   39 :         int T_rows;
##   40 :         int N;
##   41 :         double nu;
##   42 :         double tau;
##   43 :         vector_d eta;
##   44 :         matrix_d Ret;
##   45 :         matrix_d omega;
##   46 : public:
##   47 :     model99440854acd_2af95bb1eb403f4e3d61982a99281fba(rstan::io::rlist_ref_var_context& context__,
##   48 :         std::ostream* pstream__ = 0)
##   49 :         : model_base_crtp(0) {
##   50 :         ctor_body(context__, 0, pstream__);
##   51 :     }
##   52 : 
##   53 :     model99440854acd_2af95bb1eb403f4e3d61982a99281fba(stan::io::var_context& context__,
##   54 :         unsigned int random_seed__,
##   55 :         std::ostream* pstream__ = 0)
##   56 :         : model_base_crtp(0) {
##   57 :         ctor_body(context__, random_seed__, pstream__);
##   58 :     }
##   59 : 
##   60 :     void ctor_body(stan::io::var_context& context__,
##   61 :                    unsigned int random_seed__,
##   62 :                    std::ostream* pstream__) {
##   63 :         typedef double local_scalar_t__;
##   64 : 
##   65 :         boost::ecuyer1988 base_rng__ =
##   66 :           stan::services::util::create_rng(random_seed__, 0);
##   67 :         (void) base_rng__;  // suppress unused var warning
##   68 : 
##   69 :         current_statement_begin__ = -1;
##   70 : 
##   71 :         static const char* function__ = "model99440854acd_2af95bb1eb403f4e3d61982a99281fba_namespace::model99440854acd_2af95bb1eb403f4e3d61982a99281fba";
##   72 :         (void) function__;  // dummy to suppress unused var warning
##   73 :         size_t pos__;
##   74 :         (void) pos__;  // dummy to suppress unused var warning
##   75 :         std::vector<int> vals_i__;
##   76 :         std::vector<double> vals_r__;
##   77 :         local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
##   78 :         (void) DUMMY_VAR__;  // suppress unused var warning
##   79 : 
##   80 :         try {
##   81 :             // initialize data block variables from context__
##   82 :             current_statement_begin__ = 4;
##   83 :             context__.validate_dims("data initialization", "T_rows", "int", context__.to_vec());
##   84 :             T_rows = int(0);
##   85 :             vals_i__ = context__.vals_i("T_rows");
##   86 :             pos__ = 0;
##   87 :             T_rows = vals_i__[pos__++];
##   88 :             check_greater_or_equal(function__, "T_rows", T_rows, 5);
##   89 : 
##   90 :             current_statement_begin__ = 5;
##   91 :             context__.validate_dims("data initialization", "N", "int", context__.to_vec());
##   92 :             N = int(0);
##   93 :             vals_i__ = context__.vals_i("N");
##   94 :             pos__ = 0;
##   95 :             N = vals_i__[pos__++];
##   96 :             check_greater_or_equal(function__, "N", N, 2);
##   97 : 
##   98 :             current_statement_begin__ = 6;
##   99 :             context__.validate_dims("data initialization", "nu", "double", context__.to_vec());
##  100 :             nu = double(0);
##  101 :             vals_r__ = context__.vals_r("nu");
##  102 :             pos__ = 0;
##  103 :             nu = vals_r__[pos__++];
##  104 :             check_greater_or_equal(function__, "nu", nu, (N - 1));
##  105 : 
##  106 :             current_statement_begin__ = 7;
##  107 :             context__.validate_dims("data initialization", "tau", "double", context__.to_vec());
##  108 :             tau = double(0);
##  109 :             vals_r__ = context__.vals_r("tau");
##  110 :             pos__ = 0;
##  111 :             tau = vals_r__[pos__++];
##  112 :             check_greater_or_equal(function__, "tau", tau, 0);
##  113 : 
##  114 :             current_statement_begin__ = 8;
##  115 :             validate_non_negative_index("eta", "N", N);
##  116 :             context__.validate_dims("data initialization", "eta", "vector_d", context__.to_vec(N));
##  117 :             eta = Eigen::Matrix<double, Eigen::Dynamic, 1>(N);
##  118 :             vals_r__ = context__.vals_r("eta");
##  119 :             pos__ = 0;
##  120 :             size_t eta_j_1_max__ = N;
##  121 :             for (size_t j_1__ = 0; j_1__ < eta_j_1_max__; ++j_1__) {
##  122 :                 eta(j_1__) = vals_r__[pos__++];
##  123 :             }
##  124 : 
##  125 :             current_statement_begin__ = 9;
##  126 :             validate_non_negative_index("Ret", "T_rows", T_rows);
##  127 :             validate_non_negative_index("Ret", "N", N);
##  128 :             context__.validate_dims("data initialization", "Ret", "matrix_d", context__.to_vec(T_rows,N));
##  129 :             Ret = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>(T_rows, N);
##  130 :             vals_r__ = context__.vals_r("Ret");
##  131 :             pos__ = 0;
##  132 :             size_t Ret_j_2_max__ = N;
##  133 :             size_t Ret_j_1_max__ = T_rows;
##  134 :             for (size_t j_2__ = 0; j_2__ < Ret_j_2_max__; ++j_2__) {
##  135 :                 for (size_t j_1__ = 0; j_1__ < Ret_j_1_max__; ++j_1__) {
##  136 :                     Ret(j_1__, j_2__) = vals_r__[pos__++];
##  137 :                 }
##  138 :             }
##  139 : 
##  140 :             current_statement_begin__ = 10;
##  141 :             validate_non_negative_index("omega", "N", N);
##  142 :             validate_non_negative_index("omega", "N", N);
##  143 :             context__.validate_dims("data initialization", "omega", "matrix_d", context__.to_vec(N,N));
##  144 :             omega = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>(N, N);
##  145 :             vals_r__ = context__.vals_r("omega");
##  146 :             pos__ = 0;
##  147 :             size_t omega_j_2_max__ = N;
##  148 :             size_t omega_j_1_max__ = N;
##  149 :             for (size_t j_2__ = 0; j_2__ < omega_j_2_max__; ++j_2__) {
##  150 :                 for (size_t j_1__ = 0; j_1__ < omega_j_1_max__; ++j_1__) {
##  151 :                     omega(j_1__, j_2__) = vals_r__[pos__++];
##  152 :                 }
##  153 :             }
##  154 :             stan::math::check_cov_matrix(function__, "omega", omega);
##  155 : 
##  156 : 
##  157 :             // initialize transformed data variables
##  158 :             // execute transformed data statements
##  159 : 
##  160 :             // validate transformed data
##  161 : 
##  162 :             // validate, set parameter ranges
##  163 :             num_params_r__ = 0U;
##  164 :             param_ranges_i__.clear();
##  165 :             current_statement_begin__ = 16;
##  166 :             validate_non_negative_index("mu", "N", N);
##  167 :             num_params_r__ += N;
##  168 :             current_statement_begin__ = 17;
##  169 :             validate_non_negative_index("sigma", "N", N);
##  170 :             validate_non_negative_index("sigma", "N", N);
##  171 :             num_params_r__ += (N + ((N * (N - 1)) / 2));
##  172 :         } catch (const std::exception& e) {
##  173 :             stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
##  174 :             // Next line prevents compiler griping about no return
##  175 :             throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
##  176 :         }
##  177 :     }
##  178 : 
##  179 :     ~model99440854acd_2af95bb1eb403f4e3d61982a99281fba() { }
##  180 : 
##  181 : 
##  182 :     void transform_inits(const stan::io::var_context& context__,
##  183 :                          std::vector<int>& params_i__,
##  184 :                          std::vector<double>& params_r__,
##  185 :                          std::ostream* pstream__) const {
##  186 :         typedef double local_scalar_t__;
##  187 :         stan::io::writer<double> writer__(params_r__, params_i__);
##  188 :         size_t pos__;
##  189 :         (void) pos__; // dummy call to supress warning
##  190 :         std::vector<double> vals_r__;
##  191 :         std::vector<int> vals_i__;
##  192 : 
##  193 :         current_statement_begin__ = 16;
##  194 :         if (!(context__.contains_r("mu")))
##  195 :             stan::lang::rethrow_located(std::runtime_error(std::string("Variable mu missing")), current_statement_begin__, prog_reader__());
##  196 :         vals_r__ = context__.vals_r("mu");
##  197 :         pos__ = 0U;
##  198 :         validate_non_negative_index("mu", "N", N);
##  199 :         context__.validate_dims("parameter initialization", "mu", "vector_d", context__.to_vec(N));
##  200 :         Eigen::Matrix<double, Eigen::Dynamic, 1> mu(N);
##  201 :         size_t mu_j_1_max__ = N;
##  202 :         for (size_t j_1__ = 0; j_1__ < mu_j_1_max__; ++j_1__) {
##  203 :             mu(j_1__) = vals_r__[pos__++];
##  204 :         }
##  205 :         try {
##  206 :             writer__.vector_unconstrain(mu);
##  207 :         } catch (const std::exception& e) {
##  208 :             stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable mu: ") + e.what()), current_statement_begin__, prog_reader__());
##  209 :         }
##  210 : 
##  211 :         current_statement_begin__ = 17;
##  212 :         if (!(context__.contains_r("sigma")))
##  213 :             stan::lang::rethrow_located(std::runtime_error(std::string("Variable sigma missing")), current_statement_begin__, prog_reader__());
##  214 :         vals_r__ = context__.vals_r("sigma");
##  215 :         pos__ = 0U;
##  216 :         validate_non_negative_index("sigma", "N", N);
##  217 :         validate_non_negative_index("sigma", "N", N);
##  218 :         context__.validate_dims("parameter initialization", "sigma", "matrix_d", context__.to_vec(N,N));
##  219 :         Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> sigma(N, N);
##  220 :         size_t sigma_j_2_max__ = N;
##  221 :         size_t sigma_j_1_max__ = N;
##  222 :         for (size_t j_2__ = 0; j_2__ < sigma_j_2_max__; ++j_2__) {
##  223 :             for (size_t j_1__ = 0; j_1__ < sigma_j_1_max__; ++j_1__) {
##  224 :                 sigma(j_1__, j_2__) = vals_r__[pos__++];
##  225 :             }
##  226 :         }
##  227 :         try {
##  228 :             writer__.cov_matrix_unconstrain(sigma);
##  229 :         } catch (const std::exception& e) {
##  230 :             stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable sigma: ") + e.what()), current_statement_begin__, prog_reader__());
##  231 :         }
##  232 : 
##  233 :         params_r__ = writer__.data_r();
##  234 :         params_i__ = writer__.data_i();
##  235 :     }
##  236 : 
##  237 :     void transform_inits(const stan::io::var_context& context,
##  238 :                          Eigen::Matrix<double, Eigen::Dynamic, 1>& params_r,
##  239 :                          std::ostream* pstream__) const {
##  240 :       std::vector<double> params_r_vec;
##  241 :       std::vector<int> params_i_vec;
##  242 :       transform_inits(context, params_i_vec, params_r_vec, pstream__);
##  243 :       params_r.resize(params_r_vec.size());
##  244 :       for (int i = 0; i < params_r.size(); ++i)
##  245 :         params_r(i) = params_r_vec[i];
##  246 :     }
##  247 : 
##  248 : 
##  249 :     template <bool propto__, bool jacobian__, typename T__>
##  250 :     T__ log_prob(std::vector<T__>& params_r__,
##  251 :                  std::vector<int>& params_i__,
##  252 :                  std::ostream* pstream__ = 0) const {
##  253 : 
##  254 :         typedef T__ local_scalar_t__;
##  255 : 
##  256 :         local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
##  257 :         (void) DUMMY_VAR__;  // dummy to suppress unused var warning
##  258 : 
##  259 :         T__ lp__(0.0);
##  260 :         stan::math::accumulator<T__> lp_accum__;
##  261 :         try {
##  262 :             stan::io::reader<local_scalar_t__> in__(params_r__, params_i__);
##  263 : 
##  264 :             // model parameters
##  265 :             current_statement_begin__ = 16;
##  266 :             Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> mu;
##  267 :             (void) mu;  // dummy to suppress unused var warning
##  268 :             if (jacobian__)
##  269 :                 mu = in__.vector_constrain(N, lp__);
##  270 :             else
##  271 :                 mu = in__.vector_constrain(N);
##  272 : 
##  273 :             current_statement_begin__ = 17;
##  274 :             Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, Eigen::Dynamic> sigma;
##  275 :             (void) sigma;  // dummy to suppress unused var warning
##  276 :             if (jacobian__)
##  277 :                 sigma = in__.cov_matrix_constrain(N, lp__);
##  278 :             else
##  279 :                 sigma = in__.cov_matrix_constrain(N);
##  280 : 
##  281 :             // transformed parameters
##  282 :             current_statement_begin__ = 23;
##  283 :             validate_non_negative_index("sigma_scaled", "N", N);
##  284 :             validate_non_negative_index("sigma_scaled", "N", N);
##  285 :             Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, Eigen::Dynamic> sigma_scaled(N, N);
##  286 :             stan::math::initialize(sigma_scaled, DUMMY_VAR__);
##  287 :             stan::math::fill(sigma_scaled, DUMMY_VAR__);
##  288 : 
##  289 :             // transformed parameters block statements
##  290 :             current_statement_begin__ = 24;
##  291 :             stan::math::assign(sigma_scaled, multiply((1 / tau), sigma));
##  292 : 
##  293 :             // validate transformed parameters
##  294 :             const char* function__ = "validate transformed params";
##  295 :             (void) function__;  // dummy to suppress unused var warning
##  296 : 
##  297 :             current_statement_begin__ = 23;
##  298 :             size_t sigma_scaled_j_1_max__ = N;
##  299 :             size_t sigma_scaled_j_2_max__ = N;
##  300 :             for (size_t j_1__ = 0; j_1__ < sigma_scaled_j_1_max__; ++j_1__) {
##  301 :                 for (size_t j_2__ = 0; j_2__ < sigma_scaled_j_2_max__; ++j_2__) {
##  302 :                     if (stan::math::is_uninitialized(sigma_scaled(j_1__, j_2__))) {
##  303 :                         std::stringstream msg__;
##  304 :                         msg__ << "Undefined transformed parameter: sigma_scaled" << "(" << j_1__ << ", " << j_2__ << ")";
##  305 :                         stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable sigma_scaled: ") + msg__.str()), current_statement_begin__, prog_reader__());
##  306 :                     }
##  307 :                 }
##  308 :             }
##  309 :             stan::math::check_cov_matrix(function__, "sigma_scaled", sigma_scaled);
##  310 : 
##  311 : 
##  312 :             // model body
##  313 : 
##  314 :             current_statement_begin__ = 31;
##  315 :             lp_accum__.add(inv_wishart_log(sigma, nu, omega));
##  316 :             current_statement_begin__ = 32;
##  317 :             lp_accum__.add(multi_normal_log(mu, eta, sigma_scaled));
##  318 :             current_statement_begin__ = 33;
##  319 :             for (int t = 1; t <= T_rows; ++t) {
##  320 : 
##  321 :                 current_statement_begin__ = 34;
##  322 :                 lp_accum__.add(multi_normal_log(get_base1(Ret, t, "Ret", 1), mu, sigma));
##  323 :             }
##  324 : 
##  325 :         } catch (const std::exception& e) {
##  326 :             stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
##  327 :             // Next line prevents compiler griping about no return
##  328 :             throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
##  329 :         }
##  330 : 
##  331 :         lp_accum__.add(lp__);
##  332 :         return lp_accum__.sum();
##  333 : 
##  334 :     } // log_prob()
##  335 : 
##  336 :     template <bool propto, bool jacobian, typename T_>
##  337 :     T_ log_prob(Eigen::Matrix<T_,Eigen::Dynamic,1>& params_r,
##  338 :                std::ostream* pstream = 0) const {
##  339 :       std::vector<T_> vec_params_r;
##  340 :       vec_params_r.reserve(params_r.size());
##  341 :       for (int i = 0; i < params_r.size(); ++i)
##  342 :         vec_params_r.push_back(params_r(i));
##  343 :       std::vector<int> vec_params_i;
##  344 :       return log_prob<propto,jacobian,T_>(vec_params_r, vec_params_i, pstream);
##  345 :     }
##  346 : 
##  347 : 
##  348 :     void get_param_names(std::vector<std::string>& names__) const {
##  349 :         names__.resize(0);
##  350 :         names__.push_back("mu");
##  351 :         names__.push_back("sigma");
##  352 :         names__.push_back("sigma_scaled");
##  353 :     }
##  354 : 
##  355 : 
##  356 :     void get_dims(std::vector<std::vector<size_t> >& dimss__) const {
##  357 :         dimss__.resize(0);
##  358 :         std::vector<size_t> dims__;
##  359 :         dims__.resize(0);
##  360 :         dims__.push_back(N);
##  361 :         dimss__.push_back(dims__);
##  362 :         dims__.resize(0);
##  363 :         dims__.push_back(N);
##  364 :         dims__.push_back(N);
##  365 :         dimss__.push_back(dims__);
##  366 :         dims__.resize(0);
##  367 :         dims__.push_back(N);
##  368 :         dims__.push_back(N);
##  369 :         dimss__.push_back(dims__);
##  370 :     }
##  371 : 
##  372 :     template <typename RNG>
##  373 :     void write_array(RNG& base_rng__,
##  374 :                      std::vector<double>& params_r__,
##  375 :                      std::vector<int>& params_i__,
##  376 :                      std::vector<double>& vars__,
##  377 :                      bool include_tparams__ = true,
##  378 :                      bool include_gqs__ = true,
##  379 :                      std::ostream* pstream__ = 0) const {
##  380 :         typedef double local_scalar_t__;
##  381 : 
##  382 :         vars__.resize(0);
##  383 :         stan::io::reader<local_scalar_t__> in__(params_r__, params_i__);
##  384 :         static const char* function__ = "model99440854acd_2af95bb1eb403f4e3d61982a99281fba_namespace::write_array";
##  385 :         (void) function__;  // dummy to suppress unused var warning
##  386 : 
##  387 :         // read-transform, write parameters
##  388 :         Eigen::Matrix<double, Eigen::Dynamic, 1> mu = in__.vector_constrain(N);
##  389 :         size_t mu_j_1_max__ = N;
##  390 :         for (size_t j_1__ = 0; j_1__ < mu_j_1_max__; ++j_1__) {
##  391 :             vars__.push_back(mu(j_1__));
##  392 :         }
##  393 : 
##  394 :         Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> sigma = in__.cov_matrix_constrain(N);
##  395 :         size_t sigma_j_2_max__ = N;
##  396 :         size_t sigma_j_1_max__ = N;
##  397 :         for (size_t j_2__ = 0; j_2__ < sigma_j_2_max__; ++j_2__) {
##  398 :             for (size_t j_1__ = 0; j_1__ < sigma_j_1_max__; ++j_1__) {
##  399 :                 vars__.push_back(sigma(j_1__, j_2__));
##  400 :             }
##  401 :         }
##  402 : 
##  403 :         double lp__ = 0.0;
##  404 :         (void) lp__;  // dummy to suppress unused var warning
##  405 :         stan::math::accumulator<double> lp_accum__;
##  406 : 
##  407 :         local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
##  408 :         (void) DUMMY_VAR__;  // suppress unused var warning
##  409 : 
##  410 :         if (!include_tparams__ && !include_gqs__) return;
##  411 : 
##  412 :         try {
##  413 :             // declare and define transformed parameters
##  414 :             current_statement_begin__ = 23;
##  415 :             validate_non_negative_index("sigma_scaled", "N", N);
##  416 :             validate_non_negative_index("sigma_scaled", "N", N);
##  417 :             Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> sigma_scaled(N, N);
##  418 :             stan::math::initialize(sigma_scaled, DUMMY_VAR__);
##  419 :             stan::math::fill(sigma_scaled, DUMMY_VAR__);
##  420 : 
##  421 :             // do transformed parameters statements
##  422 :             current_statement_begin__ = 24;
##  423 :             stan::math::assign(sigma_scaled, multiply((1 / tau), sigma));
##  424 : 
##  425 :             if (!include_gqs__ && !include_tparams__) return;
##  426 :             // validate transformed parameters
##  427 :             const char* function__ = "validate transformed params";
##  428 :             (void) function__;  // dummy to suppress unused var warning
##  429 : 
##  430 :             current_statement_begin__ = 23;
##  431 :             stan::math::check_cov_matrix(function__, "sigma_scaled", sigma_scaled);
##  432 : 
##  433 :             // write transformed parameters
##  434 :             if (include_tparams__) {
##  435 :                 size_t sigma_scaled_j_2_max__ = N;
##  436 :                 size_t sigma_scaled_j_1_max__ = N;
##  437 :                 for (size_t j_2__ = 0; j_2__ < sigma_scaled_j_2_max__; ++j_2__) {
##  438 :                     for (size_t j_1__ = 0; j_1__ < sigma_scaled_j_1_max__; ++j_1__) {
##  439 :                         vars__.push_back(sigma_scaled(j_1__, j_2__));
##  440 :                     }
##  441 :                 }
##  442 :             }
##  443 :             if (!include_gqs__) return;
##  444 :         } catch (const std::exception& e) {
##  445 :             stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
##  446 :             // Next line prevents compiler griping about no return
##  447 :             throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
##  448 :         }
##  449 :     }
##  450 : 
##  451 :     template <typename RNG>
##  452 :     void write_array(RNG& base_rng,
##  453 :                      Eigen::Matrix<double,Eigen::Dynamic,1>& params_r,
##  454 :                      Eigen::Matrix<double,Eigen::Dynamic,1>& vars,
##  455 :                      bool include_tparams = true,
##  456 :                      bool include_gqs = true,
##  457 :                      std::ostream* pstream = 0) const {
##  458 :       std::vector<double> params_r_vec(params_r.size());
##  459 :       for (int i = 0; i < params_r.size(); ++i)
##  460 :         params_r_vec[i] = params_r(i);
##  461 :       std::vector<double> vars_vec;
##  462 :       std::vector<int> params_i_vec;
##  463 :       write_array(base_rng, params_r_vec, params_i_vec, vars_vec, include_tparams, include_gqs, pstream);
##  464 :       vars.resize(vars_vec.size());
##  465 :       for (int i = 0; i < vars.size(); ++i)
##  466 :         vars(i) = vars_vec[i];
##  467 :     }
##  468 : 
##  469 :     std::string model_name() const {
##  470 :         return "model99440854acd_2af95bb1eb403f4e3d61982a99281fba";
##  471 :     }
##  472 : 
##  473 : 
##  474 :     void constrained_param_names(std::vector<std::string>& param_names__,
##  475 :                                  bool include_tparams__ = true,
##  476 :                                  bool include_gqs__ = true) const {
##  477 :         std::stringstream param_name_stream__;
##  478 :         size_t mu_j_1_max__ = N;
##  479 :         for (size_t j_1__ = 0; j_1__ < mu_j_1_max__; ++j_1__) {
##  480 :             param_name_stream__.str(std::string());
##  481 :             param_name_stream__ << "mu" << '.' << j_1__ + 1;
##  482 :             param_names__.push_back(param_name_stream__.str());
##  483 :         }
##  484 :         size_t sigma_j_2_max__ = N;
##  485 :         size_t sigma_j_1_max__ = N;
##  486 :         for (size_t j_2__ = 0; j_2__ < sigma_j_2_max__; ++j_2__) {
##  487 :             for (size_t j_1__ = 0; j_1__ < sigma_j_1_max__; ++j_1__) {
##  488 :                 param_name_stream__.str(std::string());
##  489 :                 param_name_stream__ << "sigma" << '.' << j_1__ + 1 << '.' << j_2__ + 1;
##  490 :                 param_names__.push_back(param_name_stream__.str());
##  491 :             }
##  492 :         }
##  493 : 
##  494 :         if (!include_gqs__ && !include_tparams__) return;
##  495 : 
##  496 :         if (include_tparams__) {
##  497 :             size_t sigma_scaled_j_2_max__ = N;
##  498 :             size_t sigma_scaled_j_1_max__ = N;
##  499 :             for (size_t j_2__ = 0; j_2__ < sigma_scaled_j_2_max__; ++j_2__) {
##  500 :                 for (size_t j_1__ = 0; j_1__ < sigma_scaled_j_1_max__; ++j_1__) {
##  501 :                     param_name_stream__.str(std::string());
##  502 :                     param_name_stream__ << "sigma_scaled" << '.' << j_1__ + 1 << '.' << j_2__ + 1;
##  503 :                     param_names__.push_back(param_name_stream__.str());
##  504 :                 }
##  505 :             }
##  506 :         }
##  507 : 
##  508 :         if (!include_gqs__) return;
##  509 :     }
##  510 : 
##  511 : 
##  512 :     void unconstrained_param_names(std::vector<std::string>& param_names__,
##  513 :                                    bool include_tparams__ = true,
##  514 :                                    bool include_gqs__ = true) const {
##  515 :         std::stringstream param_name_stream__;
##  516 :         size_t mu_j_1_max__ = N;
##  517 :         for (size_t j_1__ = 0; j_1__ < mu_j_1_max__; ++j_1__) {
##  518 :             param_name_stream__.str(std::string());
##  519 :             param_name_stream__ << "mu" << '.' << j_1__ + 1;
##  520 :             param_names__.push_back(param_name_stream__.str());
##  521 :         }
##  522 :         size_t sigma_j_1_max__ = (N + ((N * (N - 1)) / 2));
##  523 :         for (size_t j_1__ = 0; j_1__ < sigma_j_1_max__; ++j_1__) {
##  524 :             param_name_stream__.str(std::string());
##  525 :             param_name_stream__ << "sigma" << '.' << j_1__ + 1;
##  526 :             param_names__.push_back(param_name_stream__.str());
##  527 :         }
##  528 : 
##  529 :         if (!include_gqs__ && !include_tparams__) return;
##  530 : 
##  531 :         if (include_tparams__) {
##  532 :             size_t sigma_scaled_j_1_max__ = (N + ((N * (N - 1)) / 2));
##  533 :             for (size_t j_1__ = 0; j_1__ < sigma_scaled_j_1_max__; ++j_1__) {
##  534 :                 param_name_stream__.str(std::string());
##  535 :                 param_name_stream__ << "sigma_scaled" << '.' << j_1__ + 1;
##  536 :                 param_names__.push_back(param_name_stream__.str());
##  537 :             }
##  538 :         }
##  539 : 
##  540 :         if (!include_gqs__) return;
##  541 :     }
##  542 : 
##  543 : }; // model
##  544 : 
##  545 : }  // namespace
##  546 : 
##  547 : typedef model99440854acd_2af95bb1eb403f4e3d61982a99281fba_namespace::model99440854acd_2af95bb1eb403f4e3d61982a99281fba stan_model;
##  548 : 
##  549 : #ifndef USING_R
##  550 : 
##  551 : stan::model::model_base& new_model(
##  552 :         stan::io::var_context& data_context,
##  553 :         unsigned int seed,
##  554 :         std::ostream* msg_stream) {
##  555 :   stan_model* m = new stan_model(data_context, seed, msg_stream);
##  556 :   return *m;
##  557 : }
##  558 : 
##  559 : #endif
##  560 : 
##  561 : 
##  562 : 
##  563 : #include <rstan_next/stan_fit.hpp>
##  564 : 
##  565 : struct stan_model_holder {
##  566 :     stan_model_holder(rstan::io::rlist_ref_var_context rcontext,
##  567 :                       unsigned int random_seed)
##  568 :     : rcontext_(rcontext), random_seed_(random_seed)
##  569 :      {
##  570 :      }
##  571 : 
##  572 :    //stan::math::ChainableStack ad_stack;
##  573 :    rstan::io::rlist_ref_var_context rcontext_;
##  574 :    unsigned int random_seed_;
##  575 : };
##  576 : 
##  577 : Rcpp::XPtr<stan::model::model_base> model_ptr(stan_model_holder* smh) {
##  578 :   Rcpp::XPtr<stan::model::model_base> model_instance(new stan_model(smh->rcontext_, smh->random_seed_), true);
##  579 :   return model_instance;
##  580 : }
##  581 : 
##  582 : Rcpp::XPtr<rstan::stan_fit_base> fit_ptr(stan_model_holder* smh) {
##  583 :   return Rcpp::XPtr<rstan::stan_fit_base>(new rstan::stan_fit(model_ptr(smh), smh->random_seed_), true);
##  584 : }
##  585 : 
##  586 : std::string model_name(stan_model_holder* smh) {
##  587 :   return model_ptr(smh).get()->model_name();
##  588 : }
##  589 : 
##  590 : RCPP_MODULE(stan_fit4model99440854acd_2af95bb1eb403f4e3d61982a99281fba_mod){
##  591 :   Rcpp::class_<stan_model_holder>("stan_fit4model99440854acd_2af95bb1eb403f4e3d61982a99281fba")
##  592 :   .constructor<rstan::io::rlist_ref_var_context, unsigned int>()
##  593 :   .method("model_ptr", &model_ptr)
##  594 :   .method("fit_ptr", &fit_ptr)
##  595 :   .method("model_name", &model_name)
##  596 :   ;
##  597 : }
##  598 : 
##  599 : 
##  600 : // declarations
##  601 : extern "C" {
##  602 : SEXP file99471e13216( ) ;
##  603 : }
##  604 : 
##  605 : // definition
##  606 : SEXP file99471e13216() {
##  607 :  return Rcpp::wrap("2af95bb1eb403f4e3d61982a99281fba");
##  608 : }
## 
## CHECKING DATA AND PREPROCESSING FOR MODEL '2af95bb1eb403f4e3d61982a99281fba' NOW.
## 
## COMPILING MODEL '2af95bb1eb403f4e3d61982a99281fba' NOW.
## 
## STARTING SAMPLER FOR MODEL '2af95bb1eb403f4e3d61982a99281fba' NOW.

##    user  system elapsed 
##    3.18    1.06 1359.21

list_of_draws <- rstan::extract(fitting)

# Extract the posterior distribution of sigma: 

sigma_posterior <- list_of_draws$sigma

# Extract the posterior distribution of mu: 

mu_posterior <- list_of_draws$mu

# Function extracts weights by Bayesian Approach: 

weights_bayesian <- function(i) {
  
  bayesian_sigma_port <- sigma_posterior[i, , ]
  
  top_mat <- cbind(2*bayesian_sigma_port, rep(1, N))
  
  bot_vec <- c(rep(1, N), 0)
  
  Am_mat <- rbind(top_mat, bot_vec)
  
  b_vec <- c(rep(0, N), 1)
  
  zm_mat <- solve(Am_mat) %*% b_vec
  
  X_min_var_baysian <- zm_mat[1:N, 1]
  
  return(X_min_var_baysian)
}


# Calculate Bayesian weights: 

number_sigmas <- length(sigma_posterior) / (N*N)

sapply(1:number_sigmas, weights_bayesian) %>% 
  t() %>% 
  colMeans() -> bayesian_weights


# Bayesian approach for asset allocation: 

bayesian_portfolio <- function(...) {
  
  bayesian_sigma_port <- sigma_posterior[1, , ]
  
  top_mat <- cbind(2*bayesian_sigma_port, rep(1, N))
  
  bot_vec <- c(rep(1, N), 0)
  
  Am_mat <- rbind(top_mat, bot_vec)
  
  b_vec <- c(rep(0, N), 1)
  
  zm_mat <- solve(Am_mat) %*% b_vec
  
  X_min_var_baysian <- zm_mat[1:N, 1]
  
  return(X_min_var_baysian)
}

2.3 A Meta-Heuristics Approach using Particle Swarm Optimization

library(pso)

# Redefine the maximum sharpe ratio function 


sharpe_func2 <- function(x){
  weight <- numeric()
  
  for (i in 1:N)
  {
    weight[i]<- x[i]
    }
  
  # Sharpe ratio numerator (Difference between return of portfolio and return of the market)
  numerator <- numeric()
  
  for (i in 1:N) {
    numerator[i] <- weight[i]*mean_ret[i]
  }
  
  # Sharpe ratio denominator
  
  denominator<- numeric()
  for (i in 1:N) {
    denominator1 <- numeric()
    for (j in 1:N) {
      denominator1[j] <- weight[i]*weight[j]*cov_mat[i,j]
    }
    denominator[i] <-sum(denominator1)
  }

sharpe_func2<- (sum(numerator) - Rf)/sum(denominator)
sharpe_func2 <- sharpe_func2 - 1e9 * (round(sum(weight),10)-1)^2

 return(sharpe_func2)
 
}


set.seed(29)
pso_port <- psoptim(par = rep(NA,4), fn =  function(x) {sharpe_func2(x)},
        lower = rep(0,4), upper = rep(1,4), 
        control = list(maxit = 10000, s = 100, maxit.stagnate = 500, fnscale = -1))

2.4 A Comparision with Weight from Markowitz Classical model

# Defining Markowitz classical optimization problem 
markowitz_portfolio <- function(...) {
  
  top_mat <- cbind(2*cov_mat, rep(1, N))
  
  bot_vec <- c(rep(1, N), 0)
  
  Am_mat <- rbind(top_mat, bot_vec)
  
  b_vec <- c(rep(0, N), 1)
  
  zm_mat <- solve(Am_mat) %*% b_vec
  
  X_min_var_markowitz <- zm_mat[1:N, 1]
  
  return(X_min_var_markowitz)
}

# Calculating Markowitz portfolio proportion

Markowitz_weight = markowitz_portfolio()

# Calculating Bayesian weight for classical model

Bayesian_classical_weight = bayesian_portfolio()

# Calculating PSO weight

PSO_weight = pso_port$par

# Scaling weight to percentage

library(scales)

## 
## Attaching package: 'scales'

## The following object is masked from 'package:purrr':
## 
##     discard

## The following object is masked from 'package:readr':
## 
##     col_factor

data.frame(stock = unique(symbols[1:4]),
          Bayesian_weight = bayes_finance_ucb$Best_Par,
          Markowitz_weight = Markowitz_weight,
          Bayesian_classical_weight = Bayesian_classical_weight,
          Particle_swarm_weight = PSO_weight) %>% 
  mutate(Bayesian_Sharpe_weight = percent(Bayesian_weight, accuracy = 0.01), 
         Markowitz_weight = percent(Markowitz_weight, accuracy = 0.01),
         Bayesian_classical_weight = percent(Bayesian_classical_weight,accuracy = 0.01),
        Particle_swarm_weight = percent(PSO_weight, accuracy = 0.01)) %>%
  arrange_all(vars(Bayesian_classical_weight, Bayesian_Sharpe_weight, Markowitz_weight, Particle_swarm_weight), desc)%>%
  select(stock, Bayesian_Sharpe_weight,Markowitz_weight, Bayesian_classical_weight, Particle_swarm_weight)

## Warning: The `...` argument of `arrange_all()` can't contain quosures. as of dplyr 0.8.3.
## Please use a one-sided formula, a function, or a function name.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

Bayesian Optimization for Portfolio Selection

Bayesian Approach

Ho Huu Binh

1 Data Preparation

2 Bayesian Approach for Portfolio Optimization Problem

2.1 Bayesian Optimization function

2.2 Bayesian Framework for Mean-Covariance model

2.3 A Meta-Heuristics Approach using Particle Swarm Optimization

2.4 A Comparision with Weight from Markowitz Classical model