improper length of one or more arguments to merge.xts

For each of the stocks, compute the mean and volatility of the monthly (simple) returns for the overall period. Plot those numbers in a scatter diagram. Then add the efficient frontier with weights constrained to be between 0 and 1, as well as the efficient frontier with weights constrained to be between 0 and 15%. Interpret the effect of the weight constraints.
Download the monthly (adjusted) price series of those stocks and ETFs for the period 2003-12-31 till 2023-08-31, together with the price series of the S&P 500.

I keep having this error : Erreur dans error(x, ...) :
improper length of one or more arguments to merge.xts
can someone help me ?


# Loading of packages ---------------------------------------------------------
library("quantmod")
library("PerformanceAnalytics")
library("quadprog")
library("PortfolioAnalytics")


# GET THE DATA ----------------------------------------------------------------

# Set the start and end date
date.start <- as.Date("2003-12-31")
date.end <- as.Date("2023-08-31")

# Define the tickers of the stocks and ETF's
tickers <- c("CAT", "XOM", "CTRA", "NFG", "GEL",   # Stocks
             "SOXX","VSMAX", "QQQ", "XLE", "XLP")  # ETF's

# Loop over the tickers and get the monthly adjusted price
monthly.close <- c()
for (ticker in tickers) {
  cat(ticker, "\n")
  price <- suppressWarnings(getSymbols(Symbols = ticker, from = date.start, to = date.end))
  adjusted.price <- get(price)[, 6]
  monthly.close <- cbind(monthly.close, xts::to.monthly(adjusted.price, indexAt = "lastof")[, 4])
}

# Rename the columns
names(monthly.close) <- tickers

# Get the data for the S&P500
getSymbols(Symbols = "^GSPC", from = date.start, to = date.end)
sp500_adjusted <- GSPC[, 6]
sp500_monthly <- xts::to.monthly(sp500_adjusted, indexAt = "lastof")[, 4]

# Calculate the returns
returns <- Return.calculate(monthly.close)
returns <- returns[(-1),]

# Returns of SP500
sp500_returns <- Return.calculate(sp500_monthly)
sp500_returns <- sp500_returns[(-1),]


# Portfolio optimization ------------------------------------------------------

# Compute the mean and volatility of the monthly returns
vMu <- colMeans(returns)
mCov <- cov(returns)

# Get the number of tickers in the dataset
N <- ncol(mCov)

# Define the grid of target return
target.mu.grid <- seq(min(vMu),max(vMu),length.out = 100 )  

# Define the weight matrix to store all optimized weights
weight.matrix     <- matrix( NA, nrow = length(target.mu.grid), ncol = N )


# Define the constraints
l <- rep(0, N)
u <- rep(1, N)
dvec <- rep(0, N)
At <- rbind( rep(1, N) , vMu , diag(rep(1, N)) , diag(rep(-1, N)) )
Amat <- t(At)


# Run the loop
for(i  in 1:length(target.mu.grid) ){
  bound <- c(1,target.mu.grid[i], l, -u)
  out <-  try(solve.QP( Dmat = mCov, dvec = dvec, 
                        Amat= Amat, bvec=bound, meq = 2)$solution, silent = TRUE)
  if(class(out)!="try-error"){
    weight.matrix[i, ] <- out
  }
}

# Remove the NAs
weight.matrix <- na.omit(weight.matrix)


# Compute the mean and sd values of optimized weights
sd.optimal <- apply(weight.matrix, 1, FUN = function(w)  sqrt(t(w) %*% mCov %*% w))
mu.optimal <- apply(weight.matrix, 1, FUN = function(w) t(w) %*% vMu)

# Compute the vector of standard deviation of the assets
vSd <- sqrt(diag(mCov))

# Plot the assets
plot(vSd, vMu, col = "gray",
     main = "Efficient frontier of selected assets",
     xlab = "Standard deviation (monthly)",
     ylab = "Average return (monthly)", 
     xlim = c(0, max(vSd) + 0.01), 
     ylim = c(0, max(vMu) + 0.005),
     las = 1)
text(vSd, vMu, labels = colnames(returns), cex = 0.7)

# Indicate in green the efficient frontier
idx <-  mu.optimal >=  mu.optimal[which.min(sd.optimal)]
lines( sd.optimal[idx], mu.optimal[idx], col = "green", lwd = 2)


# 25% weight weight constraint------------------------------------------------

# Define the weight matrix to store all optimized weights
constr.weight.matrix <- matrix(NA, nrow = length(target.mu.grid), ncol = N)

u <- rep(0.15, N) 

# Run the loop
for(i in 1:length(target.mu.grid)){
  bound <- c(1,target.mu.grid[i], l, -u)
  out <-  try(solve.QP(Dmat = mCov, dvec = dvec, 
                       Amat= Amat, bvec=bound, meq = 2)$solution, silent=TRUE)
  if(class(out)!="try-error"){
    constr.weight.matrix[i,] <- out
  }
}

# Remove the NAs
constr.weight.matrix <- na.omit(constr.weight.matrix)

# Compute the mean and sd values of optimized weights
constr.sd.optimal <- apply(constr.weight.matrix, 1, FUN = function(w) sqrt(t(w) %*% mCov %*% w))
constr.mu.optimal <- apply(constr.weight.matrix, 1, FUN = function(w) t(w) %*% vMu)

idx <- constr.mu.optimal >=  constr.mu.optimal[which.min(constr.sd.optimal)]
lines(constr.sd.optimal[idx],  constr.mu.optimal[idx], col = "blue", lwd = 2)
legend(x = "topleft",
       legend = c("weight constrain between 0 and 1", "weight constraint between 0 and 0.15"),
       col = c("green", "blue"),
       lwd = 2)



# Analyze the portfolio weights
ten_colors <- c("red", "orange", "yellow", "green", "lightblue", "blue", "violet", "purple", "pink", "grey")

barplot(t(constr.weight.matrix), beside = FALSE,
        main = "Asset weights for efficient portfolios",
        names.arg = seq(1, nrow(constr.weight.matrix)),
        xlab = "portfolio index",
        ylab = "weights",
        col = ten_colors,
        xlim = c(0,73))
legend("topright", xpd = TRUE, cex = 0.6,
       legend = rev(tickers),
       fill = rev(ten_colors))

# Compute the min and max contribution of ETF's to the portfolios
min(rowSums(constr.weight.matrix[, 6:10]))
max(rowSums(constr.weight.matrix[, 6:10]))



# BACKTESTING------------------------------------------------------------------
# without re-estimation

# We compare minimum variance, maximum Sharpe ratio and equally weighted

w.min.var <- constr.weight.matrix[which.min(constr.sd.optimal), ]
w.eq <- rep(1/N, N)
w.max.sharpe <- constr.weight.matrix[which.max(constr.mu.optimal / constr.sd.optimal), ]


# Calculate the returns for the different portfolios
returns.min.var  <- Return.portfolio(R = returns["2014-01-01/2022-02-28"],
                                     weights = w.min.var, rebalance_on = "months")

returns.eq  <- Return.portfolio(R = returns["2014-01-01/2022-02-28"],
                                weights = w.eq, rebalance_on= "months")

returns.max.sharpe <- Return.portfolio(R = returns["2014-01-01/2022-02-28"],
                                       weights = w.max.sharpe, rebalance_on = "months")


# Merge all returns into a single xts object
preturns <- merge(returns.min.var, returns.eq, returns.max.sharpe, sp500_returns["2014-01-01/2022-02-28"])

# Set the column names
colnames(preturns) <- c("Min Variance", "Equally Weighted", "Max Sharpe ratio", "S&P500")

# Evaluating the performance in terms of mean, sd and Sharpe ratio
table.AnnualizedReturns(preturns)

# Plot the performance over rolling time windows, here 24 months
charts.RollingPerformance(preturns, width = 24, legend.loc = "topleft", scale = 12)

# Plot the cumulative performance
chart.CumReturns(preturns, wealth.index = TRUE, legend.loc = "topleft",
                 main = "Cumulative performance of rebalanced portfolios") 

# Plot the annualized volatility of the portfolio
chart.RollingPerformance(R = preturns["2014-01-01/2022-02-28"],
                         width = 22, FUN = "sd.annualized", scale = 252,
                         main = "Rolling 1 month Annualized Volatility of the portfolios",
                         ylab = "", legend.loc = "topleft")


# BACKTESTING 
# with re-estimation

estim.length <- 60 # Number of months in the estimation sample
tmp <- vector(mode = "numeric", length = nrow(returns) - estim.length)
oos.perf.min.var <- oos.perf.eq <- oos.perf.max.sharpe <- tmp
n.obs <- nrow(returns)

for (j in 1:(n.obs - estim.length)){
  cat("Backtest ", j, "\n")
  
  estim.ret <- returns[j:(j + estim.length - 1), ]
  oos.ret <- returns[(j + estim.length), ]
  
  # Estimate mean and covariance matrix
  vMu <- colMeans(estim.ret)
  mCov <- cov(estim.ret)
  
  # Cardinality
  N <- ncol(mCov)
  
  # Define the grid of target return
  target.mu.grid <- seq(min(vMu),max(vMu), length.out = 100)  
  
  # Define the weight matrix to store all optimized weights
  weight.matrix     <- matrix( NA, nrow = length(target.mu.grid), ncol = N)
  
  # Define the constraints
  l <- rep(0, N)
  u <- rep(1, N)
  dvec <- rep(0, N)
  At <- rbind( rep(1, N) , vMu , diag(rep(1, N)) , diag(rep(-1, N)))
  Amat <- t(At)
  
  # Run the loop
  for(i in 1:length(target.mu.grid)){
    bound <- c(1,target.mu.grid[i] , l, -u)
    out <-  try(solve.QP( Dmat = mCov, dvec = dvec, 
                          Amat = Amat, bvec = bound, meq = 2)$solution, silent=TRUE)
    if(class(out)!="try-error"){
      weight.matrix[i,] <- out
    }
  }
  weight.matrix <- na.omit(weight.matrix)
  # Compute the mean and sd values of optimized weights
  sd.optimal <- apply(weight.matrix, 1, FUN = function(w) sqrt(t(w) %*% mCov %*% w))
  mu.optimal <- apply(weight.matrix ,1, FUN = function(w) t(w) %*% vMu)
  
  weight.min.var    <- weight.matrix[which.min(sd.optimal), ]
  # Assume risk free rate is 0
  weight.max.sharpe <- weight.matrix[which.max(mu.optimal / sd.optimal), ]
  
  oos.perf.min.var[j]    <- sum(oos.ret * weight.min.var)
  oos.perf.max.sharpe[j] <- sum(oos.ret * weight.max.sharpe)
  oos.perf.eq[j]         <- sum(oos.ret * 1 / N)
}

# Out of sample evaluation
oos.ret <- cbind(oos.perf.min.var,oos.perf.eq ,oos.perf.max.sharpe, sp500_returns["2003-12-31/2023-08-31"])
oos.ret <- xts(oos.ret, order.by = index(returns)[(estim.length + 1):nrow(returns)])
names(oos.ret)[names(oos.ret) == "sp500_adjusted.Close"] <- "sp500"


table.AnnualizedReturns(oos.ret)
# Plot the performance over rolling time windows, here 24 months
charts.RollingPerformance(R = oos.ret, width = 24, legend.loc = "topleft", scale = 12,
                          main = "24 month rolling performance of re-estimated and reballanced portfolios")


chart.CumReturns(oos.ret["2014-01-01/2023-08-31"], wealth.index = TRUE, legend.loc = "topleft",
                 main = "Cumulative return of  re-estimated and rebalanced portfolios")


# Relative performance: take 1 as reference case and plot the others relatively to them
chart.RelativePerformance(Ra = oos.ret[, 2:3]["2014-01-01/2023-08-31"],
                          Rb = oos.ret[, 1]["2014-01-01/2023-08-31"], 
                          legend.loc="topleft",
                          cex.axis = 1.3,
                          main = "Relative Performance vs minimum variance portfolio")

# Drawdown plot
chart.Drawdown(oos.ret, legend.loc = "bottomright")

# Create drowdown tables
table.Drawdowns(oos.ret[, "oos.perf.min.var"])
table.Drawdowns(oos.ret[, "oos.perf.eq"])
table.Drawdowns(oos.ret[, "oos.perf.max.sharpe"])


# Plot the annualized volatility of the portfolio
chart.RollingPerformance(R = oos.ret["2014-01-01/2023-08-31"],
                         width = 22, FUN = "sd.annualized", scale = 252,
                         legend.loc = "topleft", ylab = "",
                         main = "Rolling 1 month Annualized Volatility of the portfolios")

It works mechanically but

  1. the reprex is too long
  2. I don't know your intent
estim.length <- 60
# error comes here
# oos.ret <- xts(oos.ret, order.by = index(returns)[(estim.length + 1):nrow(returns)])
#Error in xts(oos.ret, order.by = index(returns)[(estim.length + 1):nrow(returns)]) : 
#  NROW(x) must match length(order.by)

length(oos.ret) 
length(index(returns)[(estim.length + 1)])

# Because 10 ≠ 1, so, what should order.by be?
# it and oos.ret both need to be time series and 
# a vector of equal length

# oos.ret is dim 1,10 which is one date with 10 tickers
# returns is dim 236,10, which is 236 dates with 10 tickers
# index(returns) is a vector of 236 dates only 
estim.length = 60
index(returns)[(estim.length + 1):nrow(returns)] |> length()
# to make work mechanically, this needs to be length 10,
# so let's make the argument fit
estim.length <- 226

oos.ret <- xts(oos.ret, order.by = index(returns)[estim.length])
oos.ret

I'm embarking on a project to devise and assess a specialized investment strategy in the U.S. stock market. The project's objective is to construct a balanced and efficient portfolio, analyze its performance, and offer informed investment advice. Here's an overview of my plan:

  1. Building the Portfolio: I will create a diverse portfolio consisting of 5 stocks and 5 ETFs. These securities have been chosen with specific criteria in mind: each has been listed on the public market since at least January 2004. Among the stocks, one is a company I've previously analyzed in-depth, and the others are from the same sector, ensuring sector-specific representation and coherence.
  2. Gathering Financial Data: I will acquire the monthly adjusted price data for these selected stocks and ETFs, spanning from the end of 2003 to late 2023. This comprehensive dataset will also include the price series of the S&P 500 during the same period, serving as a benchmark for comparison.
  3. Analyzing Stock Performance: For each stock, I'll calculate the average monthly returns and their volatility over the entire data range. These figures will be visually represented in a scatter plot. Additionally, I'll plot the efficient frontier – a curve representing the optimal portfolio composition for the lowest risk – under two different sets of constraints on the weights of individual assets: one allowing weights between 0 to 1, and another restricting them to a maximum of 15%. This will help me understand how these constraints affect the portfolio's risk and return profile.
  4. Optimizing the Portfolio: The focus will then shift to formulating portfolios that either minimize variance or maximize the Sharpe ratio, adhering to certain rules: all positions are long-only, the portfolio is always fully invested without any cash reserves, and no single asset exceeds 15% of the total portfolio. The aim here is to see how these conditions influence the portfolio's diversification and whether the resultant asset weights meet our expectations.
  5. Backtesting Strategies: I plan to perform an out-of-sample backtest on both the minimum variance and maximum Sharpe ratio strategies. This test will cover a period from early 2011 to late 2023, with monthly portfolio rebalancing to maintain optimized weights, and assumes no transaction costs or changes in the weighting strategy over time. The performance of these strategies will be compared against the S&P 500 and an equal-weight strategy. I'll use various charts and tables to analyze and interpret these results.
  6. Providing Investment Advice: Finally, based on the insights gained from these analyses, I will offer strategic investment advice. This advice will take into account the prevailing market conditions – such as volatility, economic growth trends, and interest rates – and how they might influence the portfolio's performance. I'll also compare the active selection strategy with the performance of the S&P 500 to determine its viability. The advice will be tailored to a specific investment horizon, ensuring it's relevant and actionable for current market scenarios.

Thank. I should have been more specific. xts returns an xts object, which is a type of time series object. What is it supposed to show in terms of the range of dates and tickers?

The xts object will cover the period from December 31, 2003, to August 31, 2023. Each row in this object corresponds to a specific point in time within this range. If the data is monthly, then each row represents the end of a month.

Tickers: The columns in the xts object represent the different tickers I analyzing. In my case, these include the tickers "CAT", "BA", "UNP", "HON", "GE" for stocks, and "SOXX", "VSMAX", "QQQ", "XLE", "XLP" for ETFs. Each column holds the time series data for one of these tickers, such as their monthly adjusted prices or calculated returns.

I'm still unclear. The two arguments to the xts() are both themselves xtx() objects ending on 2023-08-31 with the same variables having the same values:

> oos.ret
                  CAT        XOM      CTRA         NFG GEL
2023-08-31 0.06471305 0.04241818 0.0292311 0.001317979   0
                  SOXX       VSMAX         QQQ        XLE
2023-08-31 -0.05256868 -0.03718888 -0.01777525 0.01577502
                  XLP
2023-08-31 -0.0345875
> returns[236,]
                  CAT        XOM      CTRA         NFG GEL
2023-08-31 0.06471305 0.04241818 0.0292311 0.001317979   0
                  SOXX       VSMAX         QQQ        XLE
2023-08-31 -0.05256868 -0.03718888 -0.01777525 0.01577502
                  XLP
2023-08-31 -0.0345875

The return value of xts() is another xts object, but since the only row of oos.ret is the same as the last row of return, I can't understand what is the purpose of calling xts(). I understand that the same ticker values are needed, but I don't understand which rows you want to end up with.

So i suceed to correc the error but i have another one :sweat_smile:

I got this type of error : Erreur dans [.xts(oos.ret, , "oos.perf.min.var") : subscript out of bounds

Can you think someone can help me ?

# Loading of packages ---------------------------------------------------------
library("quantmod")
library("PerformanceAnalytics")
library("quadprog")
library("PortfolioAnalytics")

# GET THE DATA ----------------------------------------------------------------

# Set the start and end date
date.start <- as.Date("2003-12-31")
date.end <- as.Date("2023-08-31")

# Define the tickers of the stocks and ETF's
tickers <- c("CAT", "XOM", "CTRA", "NFG", "GEL",   # Stocks
             "SOXX", "VSMAX", "QQQ", "XLE", "XLP")  # ETF's

# Loop over the tickers and get the monthly adjusted price
monthly.close <- c()
for (ticker in tickers) {
  cat(ticker, "\n")
  price <- suppressWarnings(getSymbols(Symbols = ticker, from = date.start, to = date.end))
  adjusted.price <- get(price)[, 6]
  monthly.close <- cbind(monthly.close, xts::to.monthly(adjusted.price, indexAt = "lastof")[, 4])
}

# Rename the columns
names(monthly.close) <- tickers

# Get the data for the S&P500
getSymbols(Symbols = "^GSPC", from = date.start, to = date.end)
sp500_adjusted <- GSPC[, 6]
sp500_monthly <- xts::to.monthly(sp500_adjusted, indexAt = "lastof")[, 4]

# Calculate the returns
returns <- Return.calculate(monthly.close)
returns <- returns[(-1),]

# Returns of SP500
sp500_returns <- Return.calculate(sp500_monthly)
sp500_returns <- sp500_returns[(-1),]

# Portfolio optimization ------------------------------------------------------

# Compute the mean and volatility of the monthly returns
vMu <- colMeans(returns)
mCov <- cov(returns)

# Get the number of tickers in the dataset
N <- ncol(mCov)

# Define the grid of target return
target.mu.grid <- seq(min(vMu), max(vMu), length.out = 100)  

# Define the weight matrix to store all optimized weights
weight.matrix <- matrix(NA, nrow = length(target.mu.grid), ncol = N)

# Define the constraints
l <- rep(0, N)
u <- rep(1, N)
dvec <- rep(0, N)
At <- rbind(rep(1, N), vMu, diag(rep(1, N)), diag(rep(-1, N)))
Amat <- t(At)

# Run the loop
for (i in 1:length(target.mu.grid)) {
  bound <- c(1, target.mu.grid[i], l, -u)
  out <- try(solve.QP(Dmat = mCov, dvec = dvec, 
                      Amat = Amat, bvec = bound, meq = 2)$solution, silent = TRUE)
  if (class(out) != "try-error") {
    weight.matrix[i, ] <- out
  }
}

# Remove the NAs
weight.matrix <- na.omit(weight.matrix)

# Compute the mean and sd values of optimized weights
sd.optimal <- apply(weight.matrix, 1, FUN = function(w) sqrt(t(w) %*% mCov %*% w))
mu.optimal <- apply(weight.matrix, 1, FUN = function(w) t(w) %*% vMu)

# Compute the vector of standard deviation of the assets
vSd <- sqrt(diag(mCov))

# Plot the assets
plot(vSd, vMu, col = "gray",
     main = "Efficient frontier of selected assets",
     xlab = "Standard deviation (monthly)",
     ylab = "Average return (monthly)", 
     xlim = c(0, max(vSd) + 0.01), 
     ylim = c(0, max(vMu) + 0.005),
     las = 1)
text(vSd, vMu, labels = colnames(returns), cex = 0.7)

# Indicate in green the efficient frontier
idx <- mu.optimal >= mu.optimal[which.min(sd.optimal)]
lines(sd.optimal[idx], mu.optimal[idx], col = "green", lwd = 2)

# 25% weight constraint --------------------------------------------------------

# Define the weight matrix to store all optimized weights
constr.weight.matrix <- matrix(NA, nrow = length(target.mu.grid), ncol = N)

u <- rep(0.15, N) 

# Run the loop
for (i in 1:length(target.mu.grid)) {
  bound <- c(1, target.mu.grid[i], l, -u)
  out <- try(solve.QP(Dmat = mCov, dvec = dvec, 
                      Amat = Amat, bvec = bound, meq = 2)$solution, silent = TRUE)
  if (class(out) != "try-error") {
    constr.weight.matrix[i, ] <- out
  }
}

# Remove the NAs
constr.weight.matrix <- na.omit(constr.weight.matrix)

# Compute the mean and sd values of optimized weights
constr.sd.optimal <- apply(constr.weight.matrix, 1, FUN = function(w) sqrt(t(w) %*% mCov %*% w))
constr.mu.optimal <- apply(constr.weight.matrix, 1, FUN = function(w) t(w) %*% vMu)

idx <- constr.mu.optimal >= constr.mu.optimal[which.min(constr.sd.optimal)]
lines(constr.sd.optimal[idx], constr.mu.optimal[idx], col = "blue", lwd = 2)
legend(x = "topleft",
       legend = c("weight constrain between 0 and 1", "weight constraint between 0 and 0.15"),
       col = c("green", "blue"),
       lwd = 2)

# Analyze the portfolio weights
ten_colors <- c("red", "orange", "yellow", "green", "lightblue", "blue", "violet", "purple", "pink", "grey")

barplot(t(constr.weight.matrix), beside = FALSE,
        main = "Asset weights for efficient portfolios",
        names.arg = seq(1, nrow(constr.weight.matrix)),
        xlab = "portfolio index",
        ylab = "weights",
        col = ten_colors,
        xlim = c(0, 73))
legend("topright", xpd = TRUE, cex = 0.6,
       legend = rev(tickers),
       fill = rev(ten_colors))

# Compute the min and max contribution of ETF's to the portfolios
min(rowSums(constr.weight.matrix[, 6:10]))
max(rowSums(constr.weight.matrix[, 6:10]))

# BACKTESTING ------------------------------------------------------------------
# without re-estimation

# We compare minimum variance, maximum Sharpe ratio and equally weighted

w.min.var <- constr.weight.matrix[which.min(constr.sd.optimal), ]
w.eq <- rep(1/N, N)
w.max.sharpe <- constr.weight.matrix[which.max(constr.mu.optimal / constr.sd.optimal), ]

# Calculate the returns for the different portfolios
returns.min.var  <- Return.portfolio(R = returns["2011-01-01/2023-08-31"],
                                     weights = w.min.var, rebalance_on = "months")

returns.eq  <- Return.portfolio(R = returns["2011-01-01/2023-08-31"],
                                weights = w.eq, rebalance_on = "months")

returns.max.sharpe <- Return.portfolio(R = returns["2011-01-01/2023-08-31"],
                                       weights = w.max.sharpe, rebalance_on = "months")

# Merge all returns into a single xts object
preturns <- merge(returns.min.var, returns.eq, returns.max.sharpe, sp500_returns["2011-01-01/2023-08-31"])

# Set the column names
colnames(preturns) <- c("Min Variance", "Equally Weighted", "Max Sharpe ratio", "S&P500")

# Evaluating the performance in terms of mean, sd, and Sharpe ratio
table.AnnualizedReturns(preturns)

# Plot the performance over rolling time windows, here 24 months
charts.RollingPerformance(preturns, width = 24, legend.loc = "topleft", scale = 12)

# Plot the cumulative performance
chart.CumReturns(preturns, wealth.index = TRUE, legend.loc = "topleft",
                 main = "Cumulative performance of rebalanced portfolios") 

# Plot the annualized volatility of the portfolio
chart.RollingPerformance(R = preturns["2011-01-01/2023-08-31"],
                         width = 22, FUN = "sd.annualized", scale = 252,
                         main = "Rolling 1 month Annualized Volatility of the portfolios",
                         ylab = "", legend.loc = "topleft")

# BACKTESTING
# with re-estimation

estim.length <- 60 # Number of months in the estimation sample
tmp <- vector(mode = "numeric", length = nrow(returns) - estim.length)
oos.perf.min.var <- oos.perf.eq <- oos.perf.max.sharpe <- tmp
n.obs <- nrow(returns)

for (j in 1:(n.obs - estim.length)){
  cat("Backtest ", j, "\n")
  
  estim.ret <- returns[j:(j + estim.length - 1), ]
  oos.ret <- returns[(j + estim.length), ]
  
  # Estimate mean and covariance matrix
  vMu <- colMeans(estim.ret)
  mCov <- cov(estim.ret)
  
  # Cardinality
  N <- ncol(mCov)
  
  # Define the grid of target return
  target.mu.grid <- seq(min(vMu), max(vMu), length.out = 100)  
  
  # Define the weight matrix to store all optimized weights
  weight.matrix <- matrix(NA, nrow = length(target.mu.grid), ncol = N)
  
  # Define the constraints
  l <- rep(0, N)
  u <- rep(1, N)
  dvec <- rep(0, N)
  At <- rbind(rep(1, N), vMu, diag(rep(1, N)), diag(rep(-1, N)))
  Amat <- t(At)
  
  # Run the loop
  for (i in 1:length(target.mu.grid)) {
    bound <- c(1, target.mu.grid[i], l, -u)
    out <- try(solve.QP(Dmat = mCov, dvec = dvec, 
                        Amat = Amat, bvec = bound, meq = 2)$solution, silent = TRUE)
    if (class(out) != "try-error") {
      weight.matrix[i, ] <- out
    }
  }
  weight.matrix <- na.omit(weight.matrix)
  # Compute the mean and sd values of optimized weights
  sd.optimal <- apply(weight.matrix, 1, FUN = function(w) sqrt(t(w) %*% mCov %*% w))
  mu.optimal <- apply(weight.matrix, 1, FUN = function(w) t(w) %*% vMu)
  
  weight.min.var    <- weight.matrix[which.min(sd.optimal), ]
  # Assume the risk-free rate is 0
  weight.max.sharpe <- weight.matrix[which.max(mu.optimal / sd.optimal), ]
  
  oos.perf.min.var[j]    <- sum(oos.ret * weight.min.var)
  oos.perf.max.sharpe[j] <- sum(oos.ret * weight.max.sharpe)
  oos.perf.eq[j]         <- sum(oos.ret * 1 / N)
}


# Out of sample evaluation
oos.ret <- cbind(oos.perf.min.var, oos.perf.eq, oos.perf.max.sharpe)
oos.ret <- merge(oos.ret, sp500_returns["2014-01-01/2023-08-31"])

# Ensure the lengths match
oos_dates <- index(returns)[(estim.length + 1):n.obs]
oos.ret <- merge(xts(oos.perf.min.var, order.by = oos_dates),
                 xts(oos.perf.eq, order.by = oos_dates),
                 xts(oos.perf.max.sharpe, order.by = oos_dates),
                 xts(sp500_returns["2014-01-01/2023-08-31"], order.by = index(sp500_returns["2014-01-01/2023-08-31"])),
                 all = FALSE)

# Filter the relevant time period
oos.ret <- window(oos.ret, start = as.Date("2014-01-01"), end = as.Date("2023-08-31"))

# Plot the performance over rolling time windows, here 24 months
charts.RollingPerformance(R = oos.ret, width = 24, legend.loc = "topleft", scale = 12,
                          main = "24-month rolling performance of re-estimated and rebalanced portfolios")

chart.CumReturns(oos.ret["2014-01-01/2023-08-31"], wealth.index = TRUE, legend.loc = "topleft",
                 main = "Cumulative return of re-estimated and rebalanced portfolios")

# Relative performance: take 1 as reference case and plot the others relatively to them
chart.RelativePerformance(Ra = oos.ret[, 2:3]["2014-01-01/2023-08-31"],
                          Rb = oos.ret[, 1]["2014-01-01/2023-08-31"], 
                          legend.loc = "topleft",
                          cex.axis = 1.3,
                          main = "Relative Performance vs minimum variance portfolio")

# Drawdown plot
chart.Drawdown(oos.ret, legend.loc = "bottomright")

# Create drawdown tables
table.Drawdowns(oos.ret[, "oos.perf.min.var"])
table.Drawdowns(oos.ret[, "oos.perf.eq"])
table.Drawdowns(oos.ret[, "oos.perf.max.sharpe"])


# Plot the annualized volatility of the portfolio
chart.RollingPerformance(R = oos.ret["2014-01-01/2023-08-31"],
                         width = 22, FUN = "sd.annualized", scale = 252,
                         legend.loc = "topleft", ylab = "",
                         main = "Rolling 1 month Annualized Volatility of the portfolios")