Quantitative Analytics: Optimal Portfolio Allocation

juliomaranho Finance, Optimal Portfolio Allocation, Portfolio March 19, 2020March 19, 2020 13 Minutes

From: https://lf0.com/post/portfolio-optimisation-model/portfolio-optimisation-r/

Introduction:

The literature in portfolio optimisation has been around for decades. In this post I cover a number of traditional portfolio optimisation models. The general aim is to select a portfolio of assets out of a set of all possible portfolios being considered with a defined objective function.

The data:

The data is collected using the tidyquant() package’s tq_get() function. I then convert the daily asset prices to daily log returns using the periodReturn function from the quantmod() package. Next I construct lists of 6 months worth of daily returns using the rolling_origin() function from the rsample() package. The objective is to compute on a rolling basis the 6 month mean returns mus and the 6 month covariance matrices Sigmas on the training sets (i.e. 6 months) and apply them on the test sets (i.e. 1 month later) – monthly rebalancing.

Just as with the returns data, the same is applied to the monthly prices data.

start_date <- "2013-01-01"
end_date <- "2017-08-31"
symbols <- c("AAPL", "ABBV", "A",  "APD", "AA", "CF", "NVDA", "HOG", "WMT", "AMZN"
            #,"MSFT", "F", "INTC", "ADBE", "AMG", "AKAM", "ALB", "ALK"
            )

portfolio_prices <- tq_get(
  symbols,
  from = start_date,
  to = end_date,
) %>% 
  select(symbol, date, adjusted)


portfolio_monthly_returns <- portfolio_prices %>% 
  group_by(symbol) %>% 
  tq_transmute(
    select = adjusted,
    mutate_fun = periodReturn,
    period = "daily",                              
    type = "log",
  ) %>% 
  pivot_wider(names_from = symbol, values_from = daily.returns) %>% 
  mutate(year_month = yearmonth(date)) %>% 
  nest(-year_month) %>% 
  rolling_origin(
    initial = 6,
    assess = 1,
    skip = 0,
    cumulative = FALSE)


portfolio_monthly_prices <- portfolio_prices %>%
  pivot_wider(names_from = symbol, values_from = adjusted) %>% 
  mutate(year_month = yearmonth(date)) %>% 
  nest(-year_month) %>% 
  rolling_origin(
    initial = 6,
    assess = 1,
    skip = 0,
    cumulative = FALSE)

ListNamesDates <- map(portfolio_monthly_returns$splits, ~assessment(.x) %>%
                        select(year_month)) %>%
  plyr::ldply(., data.frame) %>% 
  pull(year_month)

#########

# Goal is to define the mean and covariance matrix on the training sets and apply them on the test sets - monthly rebalancing

mus <- map(portfolio_monthly_returns$splits, ~ analysis(.x) %>% 
             unnest(data) %>%
             select(-year_month, -date) %>%
             colMeans) %>%
  setNames(c(ListNamesDates))

Sigmas <- map(portfolio_monthly_returns$splits, ~ analysis(.x) %>%
                unnest() %>%
                tk_xts(., date_var = date) %>%
                cov(.)) %>% 
  setNames(c(ListNamesDates))

Price and Returns data

The returns data for the first split looks like the following:

AAPL	ABBV	A	APD	AA	CF	NVDA	HOG	WMT	AMZN
0.000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.000000000	0.0000000000	0.000000000	0.0000000000	0.0000000000
-0.012702708	-0.008291331	0.003575412	-0.0035002417	0.008859253	-0.004740037	0.0007860485	-0.020819530	-0.0063746748	0.0045367882
-0.028249939	-0.012713395	0.019555411	0.0133515376	0.020731797	0.022152033	0.0324601827	-0.005529864	0.0037720098	0.0025886574
-0.005899428	0.002033327	-0.007259372	-0.0009230752	-0.017429792	-0.003753345	-0.0293228286	0.006958715	-0.0096029606	0.0352948721
0.002687379	-0.022004791	-0.008022622	0.0018451781	0.000000000	-0.014769161	-0.0221704307	-0.003063936	0.0027737202	-0.0077780140
-0.015752246	0.005620792	0.026649604	0.0133906716	-0.002200109	0.034376623	-0.0226730363	0.038724890	-0.0002916231	-0.0001126237

The statistics

The Mus (mean returns) data looks like:

	x
AAPL	-0.0025241
ABBV	0.0014847
A	0.0002314
APD	0.0006546
AA	-0.0010700
CF	-0.0013681
NVDA	0.0008864
HOG	0.0008061
WMT	0.0006895
AMZN	0.0006147

The Sigmas (covariance matrix) data looks like:

	AAPL	ABBV	A	APD	AA	CF	NVDA	HOG	WMT	AMZN
AAPL	0.0004166	0.0000220	0.0000112	0.0000388	0.0000528	0.0000277	0.0000408	0.0000116	-0.0000038	-0.0000366
ABBV	0.0000220	0.0003128	0.0000490	0.0000393	0.0000137	0.0000472	0.0000511	0.0000732	0.0000258	0.0000295
A	0.0000112	0.0000490	0.0002061	0.0000720	0.0000741	0.0000885	0.0000930	0.0001175	0.0000371	0.0001139
APD	0.0000388	0.0000393	0.0000720	0.0000993	0.0000523	0.0000687	0.0000578	0.0000744	0.0000107	0.0000511
AA	0.0000528	0.0000137	0.0000741	0.0000523	0.0001485	0.0000874	0.0000685	0.0000674	-0.0000012	0.0000383
CF	0.0000277	0.0000472	0.0000885	0.0000687	0.0000874	0.0002271	0.0000706	0.0000900	-0.0000112	0.0000606
NVDA	0.0000408	0.0000511	0.0000930	0.0000578	0.0000685	0.0000706	0.0002092	0.0000706	0.0000127	0.0000853
HOG	0.0000116	0.0000732	0.0001175	0.0000744	0.0000674	0.0000900	0.0000706	0.0002393	0.0000214	0.0000895
WMT	-0.0000038	0.0000258	0.0000371	0.0000107	-0.0000012	-0.0000112	0.0000127	0.0000214	0.0000682	0.0000236
AMZN	-0.0000366	0.0000295	0.0001139	0.0000511	0.0000383	0.0000606	0.0000853	0.0000895	0.0000236	0.0003118

Comparing Portfolio Optimisation

Global Minimum Variance Portfolio

The global minimum-variance portfolio $w^{gmv}$ is a portfolio of assets with gives us the lowest possible return variance or portfolio volatility. Volatility here is used as a replacement for risk, thus with less variance in volatility correlates to less risk in an asset. The portfolio focuses only on risk and ignores expected returns.

The objective function is;

Since $\Sigma$ is unknown we can estimate it as $\hat{\Sigma}$ with the covariance matrix. In which the convex solution becomes:

The objective is that we want to find the optimial weights from the model such that our risk is minimised.

The below problem consists of our $Minimisation$ problem $\min_{w} w^{T} \Sigma w$ . The quad_form function takes the quadratic form $x^T Px$ where $x$ is a vector and $P$ is a matrix or in our case ww is our weights vector and $\Sigma$ is the covariance matrix for $A_{1}, \dots, A_{5}$ . The constraints correspond to $1^Tw=1$ in which we cannot assign negative weights to our assets and that we invest all our capital in the portfolio.

We can use the Disciplined Convex Programming (CVXR) package in R, which;

Analyses the problem,
Verifies the convexity,
Converts the problem into canonical form,
Solves the problem.

We want to find the optimial weights from the model such that our risk is minimised. We can do this by solving the optimisation problem, bind the lists into a single data frame and use ggplot2 to plot the rolling one month out of sample optimal portfolio weights – based on the previous 6 months rolled mus and Sigmas.

# 1) Function: Global Minimum Variance Portfolio
GMVPportolioFunction <- function(Sigma) {
  w <- Variable(nrow(Sigma))
  problem <- Problem(               # Initialise the problem
    Minimize(                       # minimse or maximise objective function
      quad_form(
        w, Sigma                    # Model inputs, the number of weights to consider and the covariance matrix
        )
      ), 
    constraints = list(
      w >= 0,                       # First model constraint
      sum(w) == 1))                 # Second model constrain
  Solution <- solve(problem, solver="SCS")        # Solves the problem
  return(as.vector(Solution$getValue(w)))
}

# 1a) Portfolio GMVP
GMVPPortfolio <- map(
  .x = Sigmas,
  GMVPportolioFunction)

# 1b) Portfolio GMVP
GMVPPortfolioWeights <- setNames(GMVPPortfolio, ListNamesDates) %>% 
  map(.,  ~setNames(., c(symbols))) %>% 
  map_dfr(., ~bind_rows(.), .id = "date") %>% 
  mutate(date = as.Date(paste(date, "01"), format = "%Y %b %d"))

GMVPPortfolioWeights %>%
  head() %>% 
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

date	AAPL	ABBV	A	APD	AA	CF	NVDA	HOG	WMT	AMZN
2013-07-01	0.0600841	0.0335474	0.0000000	0.1487694	0.1258906	0.0727919	0.0095632	0.0000000	0.5208891	0.0284642
2013-08-01	0.1046817	0.0404000	-0.0000009	0.0917764	0.1592917	0.0299732	0.0364528	0.0000000	0.5374240	0.0000007
2013-09-01	0.1183895	0.0000005	0.0000006	0.0844453	0.1683525	0.0246108	0.0635568	0.0000002	0.5406430	0.0000004
2013-10-01	0.0975699	0.0000000	0.0024138	0.0742436	0.1166938	0.0431301	0.0968247	0.0000000	0.5691240	0.0000000
2013-11-01	0.1480922	0.0000003	0.0620406	0.0449223	-0.0000003	0.0279872	0.1613034	0.0376023	0.5014186	0.0166322
2013-12-01	0.1135468	0.0000007	0.0678773	0.0560224	-0.0000005	0.0062812	0.1195510	0.0779063	0.5588149	0.0000001

# 1c) Portfolio GMVP
GMVPPortfolioWeights %>% 
  select(date, everything()) %>% 
  pivot_longer(cols = 2:ncol(.), values_to = "weights") %>% 
  ggplot(aes(fill = name, y = weights, x = date)) + 
  geom_bar(position = "stack", stat = "identity", width = 100) +
  scale_fill_viridis(option = "magma", discrete = TRUE, name = "Asset") +
  theme_bw() +
  ggtitle("GMVP Rolling Portfolio Adjustments") +
  xlab("Date") +
  ylab("Weights")

Markowitz Portfolio

The Markowitz Mean-Variance Portfolio is constructed as follows:

We can set different risk parameters by adjusting $\lambda$ and see how the returns are affected. This can be done by running multiple optimisation problems on the data with different $\lambda$ values. Higher $\lambda$ values places emphasis on the right hand side of the equation and thus the more risk adverse the investor is.

# 2) Function: Markowitz Portfolio
MarkowitzportolioFunction <- function(mu, Sigma, lmd) {
  w <- Variable(nrow(Sigma))
  problem <- Problem(
    Maximize(t(mu) %*% w - lmd*quad_form(w, Sigma)),
    constraints = list(
      w >= 0, 
      sum(w) == 1)
  )
  Solution <- solve(problem)
  return(as.vector(Solution$getValue(w)))
}

# 2a) Portfolio Markowitz
lmd <- c(0.25, 0.5, 0.75)
MarkowitzPortfolio <- map(lmd,
                          ~map2(
                            .x = mus,
                            .y = Sigmas,
                            MarkowitzportolioFunction, lmd = .x))

# 2b) Portfolio Markowitz
MarkowitzPortfolioWeights <- setNames(MarkowitzPortfolio, lmd) %>% 
  map(., ~setNames(., c(ListNamesDates))) %>% 
  map(., ~map(., ~setNames(., c(symbols)))) %>% 
  map(., ~map_dfr(., ~bind_rows(.), .id = "date")) %>% 
  map_dfr(., ~bind_rows(.), .id = "lambda") %>% 
  mutate(date = as.Date(paste(date, "01"), format = "%Y %b %d"))

MarkowitzPortfolioWeights %>%
  head() %>% 
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

lambda	date	ABBV	A	APD	NVDA	HOG	AMZN
0.25	2013-07-01	1.0000000	0.000000	0.0000000	0.0000000	0.0000000	0
0.25	2013-08-01	0.2938168	0.000000	0.7061832	0.0000000	0.0000000	0
0.25	2013-09-01	0.0000000	0.000000	1.0000000	0.0000000	0.0000000	0
0.25	2013-10-01	0.0000000	0.000001	0.9999984	0.0000004	0.0000001	0
0.25	2013-11-01	0.0000000	0.000000	0.0000000	0.0000000	0.0000000	1
0.25	2013-12-01	0.0000000	0.000000	0.0000000	0.0000000	0.0000000	1

We can see how the risk and return is affected by the changes in the $\lambda$ values in the following plot. As the $\lambda$ values increase the less risk we take on but also we assume less returns. Low values of $\lambda$ makes us invest in a signle asset with the weight on asset $A_{5}$ , increasing the $\lambda$ values increases the weights in other assets $A_{i}$ spreading our risk.

# 2c) Portfolio Markowitz
MarkowitzPortfolioWeights %>% 
  select(date, lambda, everything()) %>% 
  pivot_longer(cols = 3:ncol(.), values_to = "weights") %>% 
  ggplot(aes(fill = name, y = weights, x = date)) + 
  geom_bar(position = "stack", stat = "identity", width = 100) +
  facet_wrap(~lambda, ncol = 1) +
  scale_fill_viridis(option = "magma", discrete = TRUE, name = "Asset") +
  theme_bw() +
  ggtitle("Markowitz Rolling Portfolio Adjustments",
          subtitle = "For different Lambda Risk Values") +
  xlab("Date") +
  ylab("Weights")

Max Sharpe Ratio

Note: I originally made this post a few months back and when I went back to it some of the parts for the Max-Sharpe Ratio portfolio did not work as before – I think this is due to the recent update in the tidyr package. I made some minor adjustments, however this part of the post is incorrect (just look at the portfolio weights plot). I hightlight in the code when I made the minor adjustments – when I have the opportunity to go back through the code more carefully I will –

# 2) Function: Max Sharpe Ratio 
#NOTE::: When we have all negative mus for a given month we obtain an NA value - the function 
# cannot find a maximum when all mus are negative compare the MaxSharpePortfolio$`2015 Sep` and
# MaxSharpePortfolio$`2015 Oct` along with the mus$`2015 Sep` and mus$`2015 Oct` we invests 99% of the portfolio
# in Sept in ABBV since its positive. In Oct ABBV mu is negative and we cannot invest at all.
# Increasing the universe of Assets may overcome this problem.
MaxSharpeportolioFunction <- function(mu, Sigma) {
  w <- Variable(nrow(Sigma))
  problem <- Problem(
    Minimize(quad_form(w, Sigma)),
    constraints = list(
      w >= 0, t(mu) %*% w == 1)
  )
  Solution <- solve(problem, solver="SCS")
  return(as.vector(Solution$getValue(w)/sum(Solution$getValue(w))))
}

# 3a) Portfolio Max Sharpe Ratio                                 
MaxSharpePortfolio <- purrr::map2(.x = mus,
                                  .y = Sigmas,
                                  .f = ~MaxSharpeportolioFunction(.x, .y))

# 3b) Portfolio Max Sharpe Ratio
  MaxSharpePortfolioWeights <- setNames(MaxSharpePortfolio, ListNamesDates) %>%
    as_tibble(.) %>%               # Made an adjustment here
    map(.,  ~setNames(., c(symbols))) %>% 
    map_dfr(., ~bind_rows(.), .id = "date") %>% 
    mutate(date = as.Date(paste(date, "01"), format = "%Y %b %d"))

MaxSharpePortfolioWeights %>%
  head() %>% 
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

date	AAPL	ABBV	A	APD	AA	CF	NVDA	HOG	WMT	AMZN
2013-07-01	0.0009088	0.1532283	0.0371040	0.2332344	-0.1010233	0.0126946	0.2083815	-0.0262255	0.5086221	-0.0269248
2013-08-01	0.0109744	0.1485344	-0.0156878	0.3292178	-0.0080488	0.0051372	0.1146053	0.0091278	0.4086911	-0.0025515
2013-09-01	0.1217206	0.1246626	0.0087004	0.3421645	0.0000000	0.0000000	0.2414533	0.0002820	0.1610166	0.0000000
2013-10-01	0.0270424	0.0000000	0.2268394	0.3427932	0.0000000	0.0000000	0.2653166	0.1283026	0.0000000	0.0097059
2013-11-01	0.2061379	0.0000000	0.1890218	0.2234739	0.0000000	0.0533254	0.0000000	0.0000000	0.0000000	0.3280410
2013-12-01	0.2252115	0.0000000	0.0306323	0.0724210	0.0000000	0.0325004	0.0000000	0.1907158	0.1192395	0.3292796

# 3c) Portfolio Max Sharpe Ratio
MaxSharpePortfolioWeights %>% 
  select(date, everything()) %>%      # Made an adjustment here
  pivot_longer(cols = 2:ncol(.), values_to = "weights") %>% 
  ggplot(aes(fill = name, y = weights, x = date)) + 
  geom_bar(position = "stack", stat = "identity", width = 100) +
  scale_fill_viridis(option = "magma", discrete = TRUE, name = "Asset") +
  theme_bw() +
  ggtitle("Maximum Sharpe Ratio Rolling Portfolio Adjustments") +
  xlab("Date") +
  ylab("Weights")

Mu Quintile Portfolio

# 1a)
mu_quintiles <- map(mus, ~ntile(., 5) %>% 
                      setNames(c(symbols))) %>% 
  map(., ~sort(., decreasing = TRUE))

# 2)       
mu_diag_sigma <- map2(
  .x = mus,
  .y = Sigmas,
  ~.x / diag(.y)) 

mu_diag_sigma_quintiles <- mu_diag_sigma %>% 
  map(., ~ntile(., 5)) %>% 
  map(., ~setNames(., c(symbols))) %>% 
  map(., ~sort(., decreasing = TRUE))

# 1b) rank the quintiles by mu
quintiles_mu <- map2(
  .x = mus,
  .y = mu_quintiles,
  ~bind_rows(.x, .y))

quintiles_mu_weights <- map(quintiles_mu, ~ .x %>%
                              mutate_all(~ case_when(. %in% c(5, 1) ~ .5, 
                                                     TRUE ~ 0)) %>%
                              pmap_dfr(., ~ 
                                         {v1 <- c(...)
                                         v2 <- if(sum(v1) > 1) replace(v1, v1 != 0, 1/sum(v1!=0)) else v1
                                         as.list(v2)}))


quintiles_mu_weights <- map(quintiles_mu_weights, ~.x[-1,]) %>% 
  map2(.x = quintiles_mu, .y = ., ~bind_rows(.x, .y))

QuintileMuWeights <- map(quintiles_mu_weights, ~.x[3, ]) %>%  # The weights are stored in row 3.
  map_dfr(., ~bind_rows(.), .id = "date") %>% 
  mutate(date = as.Date(paste(date, "01"), format = "%Y %b %d"))

QuintileMuWeights %>%
  head() %>% 
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

date	AAPL	ABBV	A	APD	AA	CF	NVDA	WMT	AMZN
2013-07-01	0.25	0.25	0.00	0.00	0.00	0.25	0.25	0.00	0.00
2013-08-01	0.00	0.25	0.00	0.25	0.25	0.25	0.00	0.00	0.00
2013-09-01	0.00	0.00	0.00	0.25	0.25	0.25	0.25	0.00	0.00
2013-10-01	0.00	0.00	0.25	0.25	0.25	0.00	0.00	0.25	0.00
2013-11-01	0.00	0.25	0.00	0.25	0.00	0.00	0.00	0.25	0.25
2013-12-01	0.25	0.00	0.00	0.00	0.00	0.00	0.25	0.25	0.25

QuintileMuWeights %>% 
    select(date, everything()) %>%
  pivot_longer(cols = 2:ncol(.), values_to = "weights") %>% 
  ggplot(aes(fill = name, y = weights, x = date)) + 
  geom_bar(position = "stack", stat = "identity", width = 100) +
  scale_fill_viridis(option = "magma", discrete = TRUE, name = "Asset") +
  theme_bw() +
  ggtitle("Quintile MU Rolling Portfolio Adjustments") +
  xlab("Date") +
  ylab("Weights")

Mu / Diag(Sigma) Portfolio

# 2a) 
mu_diag_sigma <- map2(
  .x = mus,
  .y = Sigmas,
  ~.x / diag(.y)) 

mu_diag_sigma_quintiles <- mu_diag_sigma %>% 
  map(., ~ntile(., 5)) %>% 
  map(., ~setNames(., c(symbols))) %>% 
  map(., ~sort(., decreasing = TRUE))

# 2b) rank the quintiles by mu / diag(Sigma)
quintiles_mu_diag_sigma <- map2(
  .x = mu_diag_sigma,
  .y = mu_diag_sigma_quintiles,
  ~bind_rows(.x, .y))

quintiles_mu_diag_sigma_weights <- map(quintiles_mu_diag_sigma, ~ .x %>%
                                         mutate_all(~ case_when(. %in% c(5, 1) ~ .5, 
                                                                TRUE ~ 0)) %>%
                                         pmap_dfr(., ~ 
                                                    {v1 <- c(...)
                                                    v2 <- if(sum(v1) > 1) replace(v1, v1 != 0, 1/sum(v1!=0)) else v1
                                                    as.list(v2)}))

quintiles_mu_diag_sigma_weights <- map(quintiles_mu_diag_sigma_weights, ~.x[-1,]) %>% 
  map2(.x = quintiles_mu_diag_sigma, .y = ., ~bind_rows(.x, .y))

QuintileMuDiagSigmaWeights <- map(quintiles_mu_diag_sigma_weights, ~.x[3, ]) %>% # The weights are stored in row 3.
  map_dfr(., ~bind_rows(.), .id = "date") %>% 
  mutate(date = as.Date(paste(date, "01"), format = "%Y %b %d"))

QuintileMuDiagSigmaWeights %>%
  head() %>% 
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

date	AAPL	ABBV	A	APD	AA	CF	NVDA	HOG	WMT	AMZN
2013-07-01	0.25	0.00	0.00	0.25	0.25	0.00	0.00	0.00	0.25	0.00
2013-08-01	0.00	0.00	0.00	0.25	0.25	0.25	0.00	0.00	0.25	0.00
2013-09-01	0.00	0.00	0.00	0.25	0.25	0.25	0.25	0.00	0.00	0.00
2013-10-01	0.00	0.00	0.00	0.25	0.25	0.00	0.25	0.00	0.25	0.00
2013-11-01	0.00	0.25	0.25	0.00	0.00	0.00	0.00	0.00	0.25	0.25
2013-12-01	0.00	0.00	0.00	0.00	0.25	0.00	0.25	0.25	0.25	0.00

QuintileMuDiagSigmaWeights %>% 
  select(date, everything()) %>%
  pivot_longer(cols = 2:ncol(.), values_to = "weights") %>% 
  ggplot(aes(fill = name, y = weights, x = date)) + 
  geom_bar(position = "stack", stat = "identity", width = 100) +
  scale_fill_viridis(option = "magma", discrete = TRUE, name = "Asset") +
  theme_bw() +
  ggtitle("Quintile MU Diag Sigma Rolling Portfolio Adjustments") +
  xlab("Date") +
  ylab("Weights")

Mu / Diag(sqrt(Sigma)) Portfolio

# 3a)
mu_diag_sqrt_sigma <- map2(
  .x = mus,
  .y = Sigmas,
  ~.x / sqrt(diag(.y)))  

mu_diag_sqrt_sigma_quintiles <- mu_diag_sqrt_sigma %>% 
  map(., ~ntile(., 5)) %>% 
  map(., ~setNames(., c(symbols))) %>% 
  map(., ~sort(., decreasing = TRUE))

# 3b) rank the quintiles by mu / sqrt(diag(sigma))

quintiles_mu_diag_sqrt_sigma <- map2(
  .x = mu_diag_sqrt_sigma,
  .y = mu_diag_sqrt_sigma_quintiles,
  ~bind_rows(.x, .y))

quintiles_mu_diag_sqrt_sigma_weights <- map(quintiles_mu_diag_sqrt_sigma, ~ .x %>%
                                         mutate_all(~ case_when(. %in% c(5, 1) ~ .5, 
                                                                TRUE ~ 0)) %>%
                                         pmap_dfr(., ~ 
                                                    {v1 <- c(...)
                                                    v2 <- if(sum(v1) > 1) replace(v1, v1 != 0, 1/sum(v1!=0)) else v1
                                                    as.list(v2)}))

quintiles_mu_diag_sqrt_sigma_weights <- map(quintiles_mu_diag_sqrt_sigma_weights, ~.x[-1,]) %>% 
  map2(.x = quintiles_mu_diag_sqrt_sigma, .y = ., ~bind_rows(.x, .y))

QuintileMuDiagSqrtSigmaWeights <- map(quintiles_mu_diag_sqrt_sigma_weights, ~.x[3, ]) %>% # The weights are stored in row 3.
  map_dfr(., ~bind_rows(.), .id = "date") %>% 
  mutate(date = as.Date(paste(date, "01"), format = "%Y %b %d"))

QuintileMuDiagSqrtSigmaWeights %>%
  head() %>% 
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

date	AAPL	ABBV	A	APD	AA	CF	NVDA	HOG	WMT	AMZN
2013-07-01	0.25	0.25	0.00	0.00	0.00	0.25	0.00	0.00	0.25	0.00
2013-08-01	0.00	0.00	0.00	0.25	0.25	0.25	0.00	0.00	0.25	0.00
2013-09-01	0.00	0.00	0.00	0.25	0.25	0.25	0.25	0.00	0.00	0.00
2013-10-01	0.00	0.00	0.25	0.25	0.25	0.00	0.00	0.00	0.25	0.00
2013-11-01	0.00	0.25	0.00	0.25	0.00	0.00	0.00	0.00	0.25	0.25
2013-12-01	0.00	0.00	0.00	0.00	0.25	0.00	0.25	0.25	0.00	0.25

QuintileMuDiagSqrtSigmaWeights %>% 
  select(date, everything()) %>%
  pivot_longer(cols = 2:ncol(.), values_to = "weights") %>% 
  ggplot(aes(fill = name, y = weights, x = date)) + 
  geom_bar(position = "stack", stat = "identity", width = 100) +
  scale_fill_viridis(option = "magma", discrete = TRUE, name = "Asset") +
  theme_bw() +
  ggtitle("Quintile MU Diag Sqrt(Sigma) Rolling Portfolio Adjustments") +
  xlab("Date") +
  ylab("Weights")

Analysis of Results

There are 8 different portfolio optimisation models. Putting all the data together, it’s easier to inspect the results.

############################# Assess the Returns ##################################
MonthlyReturns <- map(portfolio_monthly_prices$splits, ~assessment(.x)) %>% 
  map(., ~unnest(.x, data) %>% 
        select(-year_month) %>% 
        tk_xts(date_var = date) %>% 
        lapply(., periodReturn, period = "monthly")) %>% 
  setNames(., ListNamesDates) %>% 
  map(., ~data.frame(.)) %>% 
  bind_rows(., .id = "date") %>% 
  mutate(date = as.Date(paste(date, "01"), format = "%Y %b %d")) %>% 
  tk_xts(., date_var = date)

############# Assess the Monthly Portfolio Returns Depeninding on Weights #########

# 1d)
TS_GMVPPortfolioWeights <- GMVPPortfolioWeights %>% 
  tk_xts(., date_var = date)

PortRets_GMVP <- Return.portfolio(MonthlyReturns, weights = TS_GMVPPortfolioWeights) %>% 
  setNames(c("GMVP"))


# 2d)
TS_MarkowitzPortfolioWeights <- MarkowitzPortfolioWeights %>% 
  group_split(lambda) %>% 
  setNames(c(lmd)) %>% 
  map(., ~tk_xts(., select = -lambda, date_var = date))

MarkowitzWeightedPortfolioReturnsFunction <- function(lmd){
  WeightedReturns = Return.portfolio(MonthlyReturns, weights = TS_MarkowitzPortfolioWeights[[lmd]])
  return(fortify.zoo(WeightedReturns))
}

PortRets_Markowitz <- lapply(seq(1:length(lmd)), MarkowitzWeightedPortfolioReturnsFunction) %>% 
  setNames(c(lmd)) %>%
  bind_rows(., .id = "lmd") %>% 
  #column_to_rownames("Index") %>% 
  #select(-matches("Index")) %>% 
  #rownames_to_column("Index") %>% 
  mutate(Index = as.Date(Index),
         lmd = as.numeric(lmd)) %>% 
  pivot_wider(names_from = lmd, values_from = portfolio.returns) %>% 
  tk_xts(date_var = Index) %>% 
  setNames(c(paste("Markowitz Lambda", lmd)))

# 3d)
TS_MaxSharpePortfolioWeights <- MaxSharpePortfolioWeights %>% 
  tk_xts(., date_var = date)

PortRets_MaxSharpe <- Return.portfolio(MonthlyReturns, weights = TS_MaxSharpePortfolioWeights) %>% 
  setNames(c("MaxSharpe"))


# 4) quintile mu

TS_QuintileMuWeights <- QuintileMuWeights %>% 
  tk_xts(., date_var = date)

PortRets_Mu_Quintiles <- Return.portfolio(MonthlyReturns, weights = TS_QuintileMuWeights) %>% 
  setNames(c("Mu_Quintiles"))

# 4) quintile mu diag sigma

TS_QuintileMuDiagSigmaWeights <- QuintileMuDiagSigmaWeights %>% 
  tk_xts(., date_var = date)

PortRets_MuDiagSigma_Quintiles <- Return.portfolio(MonthlyReturns, weights = TS_QuintileMuDiagSigmaWeights) %>% 
  setNames(c("MuDiagSigma_Quintiles"))

# 4) Quintile mu diga sqrt sigma

TS_QuintileMuDiagSqrtSigmaWeights <- QuintileMuDiagSqrtSigmaWeights %>% 
  tk_xts(., date_var = date)

PortRets_MuDiagSqrtSigma_Quintiles <- Return.portfolio(MonthlyReturns, weights = TS_QuintileMuDiagSqrtSigmaWeights) %>% 
  setNames(c("MuDiagSqrtSigma_Quintiles"))

#####################################################################################

Looking at the plots the Global Minimum Variance Portfolio shows the lowest volatility in the portfolio returns. The Max Sharpe portfolio is clearly wrong here (see, previous point on Max Sharpe Ratio section). The annaulised performance metrics are also displayed at the bottom for each portfolio over the period.

########## Assess the Annualised Returns based on portfolio weights ##################
AllReturns <- cbind(PortRets_GMVP, PortRets_Markowitz, PortRets_MaxSharpe, PortRets_Mu_Quintiles,
                    PortRets_MuDiagSigma_Quintiles, PortRets_MuDiagSqrtSigma_Quintiles)
                    
library(data.table)
AllReturns %>%
  data.table(keep.rownames = TRUE) %>%
  as_tibble() %>% 
  pivot_longer(cols = 2:ncol(.)) %>% 
  ggplot(aes(x = index, y = value, color = name)) +
  geom_line() +
  facet_wrap(~name, ncol = 2) +
  ggtitle("Portfolio Returns", subtitle = "With portfolio optimised weights") +
  xlab("Returns") +
  ylab("Date") +
  theme_tq()

AnnualMetrics <- bind_cols(table.AnnualizedReturns(PortRets_GMVP) %>% 
                          rownames_to_column("Metric"),
                        table.AnnualizedReturns(PortRets_Markowitz),
                        table.AnnualizedReturns(PortRets_MaxSharpe),
                        table.AnnualizedReturns(PortRets_Mu_Quintiles),
                        table.AnnualizedReturns(PortRets_MuDiagSigma_Quintiles),
                        table.AnnualizedReturns(PortRets_MuDiagSqrtSigma_Quintiles)
                        )

AnnualMetrics %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Metric	GMVP	Markowitz Lambda 0.25	Markowitz Lambda 0.5	Markowitz Lambda 0.75	MaxSharpe	Mu_Quintiles	MuDiagSigma_Quintiles	MuDiagSqrtSigma_Quintiles
Annualized Return	0.1009	0.8441	0.8142	0.7712	0.1450	0.1939	0.2590	0.2691
Annualized Std Dev	0.1127	0.3617	0.3516	0.3353	0.1473	0.1800	0.1485	0.1535
Annualized Sharpe (Rf=0%)	0.8951	2.3340	2.3157	2.2999	0.9845	1.0774	1.7437	1.7529

Plotting the cumulative returns shows that the Markowitz lambda portfolios have the highest returns over the period, however they also have the highest standard deviation over the same period.

chart.CumReturns(AllReturns, main = "Weighted Returns by Objective Function and Risk Tolerance",
                 wealth.index = TRUE, legend.loc = "topleft", colorset = rich6equal)

charts.PerformanceSummary(AllReturns, main = "Performance for Different Weighted Assets", 
                          wealth.index = TRUE, colorset = rich6equal)

Additional

(Some additional test data I leave here for future reference)

data <- data.frame(
  A_1 = c(2.3, 0.93, 0.62, 0.74, -0.23),
  A_2 = c(0.93, 1.40, 0.2, 0.56, 0.26),
  A_3 = c(0.62, 0.2, 1.80, 0.78, -0.27),
  A_4 = c(0.74, 0.56, 0.78, 3.4, -0.56),
  A_5 = c(-0.23, 0.26, -0.27, -0.56, 2.60),
  Security = c(1, 2, 3, 4, 5),
  R_i = c(15.1, 12.5, 14.7, 9.02, 17.68)
)

Sigma <- data[, 1:5]
w <- Variable(nrow(data))
mu <- data[, 7]
n_samples = 10
lambdas <- 10^seq(-2, 3, length.out = n_samples)

ret_data <- rep(0, n_samples)
risk_data <- rep(0, n_samples)
w_data <- matrix(0, nrow = n_samples, ncol = nrow(data))

for(i in seq_along(lambdas)) {
  lambda <- lambdas[i]
  problem <- Problem(
    Maximize(
      t(mu) %*% w - lambda*quad_form(w, Sigma)
    ),
    constraints = list(w >= 0, sum(w) == 1))
  solution <- solve(problem)
  w_data[i,] <- solution$getValue(w)
  risk_data[i] <- solution$getValue(sqrt(quad_form(w, Sigma)))
  ret_data[i] <- solution$getValue(t(mu) %*% w)
}

ggplot() +
  geom_line(mapping = aes(x = risk_data, y = ret_data), color = "blue") +
  geom_point(mapping = aes(x = sqrt(diag(as.matrix(Sigma))), y = mu), color = "red")

Quantitative Analytics: Optimal Portfolio Allocation

Introduction:

The data:

Price and Returns data

The statistics

Comparing Portfolio Optimisation

Global Minimum Variance Portfolio

Markowitz Portfolio

Max Sharpe Ratio

Mu Quintile Portfolio

Mu / Diag(Sigma) Portfolio

Mu / Diag(sqrt(Sigma)) Portfolio

Analysis of Results

Additional

Published by juliomaranho

Leave a comment Cancel reply

Introduction:

The data:

Price and Returns data

The statistics

Comparing Portfolio Optimisation

Global Minimum Variance Portfolio

Markowitz Portfolio

Max Sharpe Ratio

Mu Quintile Portfolio

Mu / Diag(Sigma) Portfolio

Mu / Diag(sqrt(Sigma)) Portfolio

Analysis of Results

Additional

Share this:

Related

Published by juliomaranho

Leave a comment Cancel reply