library(tidyverse);library(lubridate);library(tbl2xts); library(glue); library(RiskPortfolios)

# Load data and look at it:

cpi <- fmxdat::cpi

Rets <- fmxdat::mvo_example_data

ref <- fmxdat::mvo_example_data_ref

constraints <- fmxdat::mvo_ex_constraints

CPI <- constraints$CPI
ZAR_dep <- constraints$ZAR_Dep

# Sample Real Returns:
Return_Forecasts <-
Rets %>% 
  gather(Alias, Returns, -date) %>% 
  left_join(., ref, by = "Alias") %>% 
  group_by(Alias) %>% 
  summarise(
                Smpl_Geom_Mean = exp(mean(log(1+Returns)))^12 - 1,
                # For your own sanity, the above is effectively the same as doing:
                # Smpl_Geom_Mean_Sanity = prod(1+Returns_Use)^(12/n())-1,
                Smpl_SD = sd(Returns) * sqrt(12)
                ) %>% 
  
  left_join(., constraints$Nominal_FC) %>% 
  
  mutate( across(c(-Alias, -Smpl_SD),  ~ifelse( !grepl("LCL", Alias), . + ZAR_dep , .))) %>% 
  mutate( across(c(-Alias, -Smpl_SD),  ~. - CPI))

Return_Forecasts <- Return_Forecasts %>% filter(Alias %in% ref$Alias)

Return_Forecasts

# If there's a missing inflation number, replace NA with:
cpireplacement <- 0.06/12

# Make all my ZAR returns real (offshore returns are already in ZAR):
Real_Rets <- 
left_join(Rets %>% gather(Alias, Returns, -date) %>% fmxdat::YM(), cpi %>% select(-date), by = "YM") %>% mutate(cpi = coalesce(cpi, cpireplacement)) %>% mutate(Returns_Use = (1+Returns)/(1+cpi)-1) %>% 
  filter(Alias %in% ref$Alias)
        
order <- ref$Alias %>% unique

# I ensure my constraints match by the order specified...
UB <- constraints$lbub %>% arrange(match(Alias, order)) %>% select(Alias, UB)
LB <- constraints$lbub %>% arrange(match(Alias, order)) %>% select(Alias, LB)

avec_lt <- 
  constraints$group_LB %>% 
  filter(Alias %in% order) %>% 
  mutate(across(-Alias, ~ifelse(is.na(.x), 0, 1))) %>% 
  # Get our ordering right:
  arrange(match(Alias, order)) %>% select(-Alias)

avec_gt <- 
  constraints$group_UB %>% 
  filter(Alias %in% order) %>% 
  mutate(across(-Alias, ~ifelse(is.na(.x), 0, 1))) %>% 
  # Get our ordering right:
  arrange(match(Alias, order)) %>% select(-Alias)

bvec_lt <- constraints$group_LB %>% select(-Alias) %>% apply(., 2, min, na.rm=T)
bvec_gt <- constraints$group_UB %>% select(-Alias) %>% apply(., 2, max, na.rm=T)

# You can now check that these constraints are valid....:
# ncol(avec_lt) == length(bvec_lt)
if( sum(avec_gt) != sum(avec_lt) ) stop("\n\nCheck const_gt and const_lt input dimensions match....If there is a value for one asset in a column, it must have an entry for both... \n\n")


mvo_constraints <- list()

mvo_constraints$order <- order
# Individual UB / LB constraints:
mvo_constraints$Indiv$UB <- UB %>% pull(UB)
mvo_constraints$Indiv$LB <- LB %>% pull(LB)

# Amat for LB and UB:::
mvo_constraints$Group$amat <- 
  bind_cols(avec_lt %>% purrr::set_names(names(avec_gt)), 
            avec_gt %>% purrr::set_names(glue("LB_{names(avec_lt)}"))) 

# bvec for LB and UB:::
mvo_constraints$Group$UB <- bvec_gt %>% purrr::set_names(glue("{names(bvec_gt)}")) 
mvo_constraints$Group$LB <- bvec_lt %>% purrr::set_names(glue("LB_{names(bvec_lt)}"))

# Long only constraint:::

mvo_constraints$Active$amat <- data.frame(long = rep(1,length(order)), short = 0, neutral = 0) %>% as.matrix()
mvo_constraints$Active$bvec <- data.frame(long = 1, short = 0, neutral = 0) %>% as.matrix()
    
# Forecast_Choice = "Default"
Forecast_Choice_List <- Return_Forecasts %>% select(-Alias, -Smpl_SD) %>% names() %>% as.list()
gamma_List <- seq(0,1,by = 0.001) %>% as.list()

LetsOptimize <- function(Real_Rets, Cov_Start_Date = ymd(19000101), Return_Forecasts, Forecast_Choice_List, mvo_constraints, gamma_List, order){
  
  # FC_Choice <- Forecast_Choice_List[[1]]
  FC_Choice <- Forecast_Choice_List

# Calculate the covariance matrix to use:
return_mat <- 
  Real_Rets %>% filter(date >= Cov_Start_Date) %>% 
  select(date, Alias, Returns_Use) %>% 
  spread(Alias, Returns_Use) %>% select(-date) %>% as.matrix()

if( any(is.na(return_mat)) ) stop("Return matrix not full... check.")

# Use Ledoit-Wolfe covariance shrinkageL
Sigma <- RiskPortfolios::covEstimation(return_mat, control = list(cov_type = "lw"))
diag(Sigma) <- diag(Sigma) + 1e-8

# Again, order matters big time here!
Sigma <- Sigma[order,order]

# Calculate forecast returns to use:

mu <- Return_Forecasts %>% select(Alias, Mu = all_of(FC_Choice)) %>% arrange(match(Alias, order)) %>% pull(Mu)
names(mu) <- order

  # Remember, this is a quadratic program solve....
safeOpt <- purrr::safely(quadprog::solve.QP)

# Loop across the different risk paramater settings:

optim_loop <- function(safeOpt, gamma_List, mvo_constraints, mu, Sigma){
  
# gamma_choice <- gamma_List[[1]]
gamma_choice <- gamma_List

# Remember::::> set UB as -1! 
groups_mat <- mvo_constraints$Group$amat %>% mutate(across(!starts_with("LB"), ~.*-1)) %>% as.matrix()

n <- length(mu)
dvec <- mu
 
# Let's run the optimizer and save:
Amat <- cbind(1, -dvec, -diag(n), diag(n), groups_mat)
bvec <- c(1, -gamma_choice, -mvo_constraints$Indiv$UB, mvo_constraints$Indiv$LB, -mvo_constraints$Group$UB, mvo_constraints$Group$LB)
opt <- safeOpt(Sigma, dvec, Amat, bvec, meq = 1)
  
opt_weights <- opt$result$solution
if(!is.null(opt_weights)) { 
  names(opt_weights) <- order
  
  return(tibble(Assets = order, Weights = round(opt_weights, 5), gamma = gamma_choice, mu = dvec, Sigma = list(Sigma)))
  
} else {
  
    return(NULL)
  
  }

  
}

Optimum_df <- gamma_List %>% map_df(~optim_loop(safeOpt, ., mvo_constraints, mu, Sigma))

if(!is.null(Optimum_df)){
  
  return(Optimum_df %>% mutate(Type = FC_Choice))
  
} else {
  
  return(NULL)
}


}

Result_all_mu <- 
Forecast_Choice_List %>% map_df(~LetsOptimize(Real_Rets, Cov_Start_Date = ymd(19000101), Return_Forecasts, ., mvo_constraints, gamma_List, order))

# So let's create and plot the efficient frontier line
TypeL <- Result_all_mu$Type %>% unique %>% as.list()
gamma_L <- Result_all_mu %>% group_split(gamma,Type)
N <- Result_all_mu$Assets %>% unique %>% length

# Also, we want to make returns nominal now again...
CPI_Add <- constraints$CPI

EF_creator <- function(gamma_L, N, CPI_Add){
  
  # gsel <- gamma_L[[3]]
  gsel <- gamma_L
  
  if(nrow(gsel) > N) stop("Multiple Optimal points...{unique(gsel$gamma)} | {unique(gsel$Type)}")
  
  Sigma <- gsel %>% slice(1) %>% select(Sigma) %>% pull(Sigma)
  
  Sigma <- Sigma[[1]][gsel$Assets,gsel$Assets]
  wt <- gsel %>% pull(Weights)
  mu <- gsel %>% mutate(mu = mu + CPI_Add) %>% pull(mu)
  
  # Now, math from above:
    tibble(
      gamma = unique(gsel$gamma),
      Eff_Mu = as.numeric(wt %*% mu),
      Eff_SD = sqrt( as.numeric( t(wt) %*% Sigma %*% wt ) ) * sqrt(12),
      Type = unique(gsel$Type)
    )
      
  }
  
  EFLine <- gamma_L %>% map_df(~EF_creator(., N, CPI_Add))