When I run a shiny app for portfolio optimization, I get this error ("constraints are inconsistent, no solution") randomly about half the time, while it seems to work correctly half the time. I'm not really sure where the code is going wrong, as I've run the specific solve.QP and quadprog related commands on the console and it works fine there. It's only when I run it on the shiny app that the issue occurs. I suspect it may have something to do with how the code is processing the input, or with how the meq is defined. But I can't tell for sure. I have looked at similar questions on stack overflow but am unsure of how exactly I should change my meq constraints or make the Dmat matrix symmetric (if those two are the issues).
Here's the code:
library(quantmod) library(lubridate) library(dplyr) library(data.table) library(quadprog) library(shiny)
Define UI for application that draws a histogram
ui <- fluidPage(
Application title
titlePanel("Robo-Advisor Shiny App"),
Sidebar with a slider inputs
fluidRow(
column(3,
numericInput(
inputId = "start",
label = "Beginning Date (yyymmdd):",
value = 20160101
),
numericInput(
inputId = "end",
label = "Ending Date (yyymmdd):",
value = 20201231
),
selectInput(
inputId = "parameter",
label = "Return optimal portfolio for given:",
choices = c("mu", "vol"),
selected = "mu"
),
numericInput(
inputId = "desired_annual_expected_return",
label = "Desired annual expected return (in decimal format):",
value = 0.2
),
numericInput(
inputId = "desired_annual_vol",
label = "Desired annual vol (in decimal format):",
value = 0.15
)
),
column(9,align="center",
fluidRow(
#splitLayout(div(plotOutput("prcPlot")), div(tableOutput("titleTable"), style = "font-size:100%"), div(tableOutput("prcTable"), style = "font-size:100%"), cellWidths = c("50%", "50%"))
div(plotOutput("prcPlot")),
div(tableOutput("titleTable")),
div(tableOutput("prcTable"))
)
)
) #close fluidRow )#close fluidPage
Define server logic required to draw a histogram
server <- function(input, output) {
#####PREPARING DATA FOR PLOT & TABLE FUNCTIONS dataPrep = reactive ({
#define variables from input boxes--------------------
startdt = ymd(input$start)
enddt = ymd(input$end)
parameter = input$parameter
d_mean = input$desired_annual_expected_return
d_sd = input$desired_annual_vol
#download and organize data---------------------------
symbolList = c("MSFT", "WMT", "AAPL", "IBM", "KO")
getSymbols(symbolList, from = startdt, to = enddt, src="yahoo") #default source is finance.yahoo.com
#Convert to dataframe
MSFT = as.data.frame(MSFT)
WMT = as.data.frame(WMT)
AAPL = as.data.frame(AAPL)
IBM = as.data.frame(IBM)
KO = as.data.frame(KO)
MSFT = to.monthly(MSFT) #converts to monthly frequency
WMT = to.monthly(WMT)
AAPL = to.monthly(AAPL)
IBM = to.monthly(IBM)
KO = to.monthly(KO)
prices = cbind(MSFT$MSFT.Adjusted, WMT$WMT.Adjusted, AAPL$AAPL.Adjusted,
IBM$IBM.Adjusted, KO$KO.Adjusted)
len = dim(prices)[1]
returns = as.data.frame(prices[2:len,] / prices[1:(len-1),]) - 1
names(returns) = c("msft", "wmt", "aapl", "ibm", "ko")
#Define the mean-variance optimization function
QPoptim = function(Eport, noshort, N, muvec, covmat){
ones = array(1,N)
Dmat = covmat
dvec = array(0,N)
Amat = cbind(muvec, ones)
b0vec = c(Eport, 1)
if(noshort==1) {
idmat = matrix(0,N,N)
diag(idmat) = 1
Amat = cbind(Amat,idmat)
b0vec = c(b0vec, array(0,N))
}
wvec = solve.QP(Dmat, dvec, Amat, b0vec, meq=2)$solution
sigport = sqrt( t(wvec) %*% covmat %*% wvec ) #returns this last value
}
#repeating the Qpoptim code partly to get only the portfolio weights
meanvarweights = function(Eport, noshort, N, muvec, covmat){
ones = array(1,N)
Dmat = covmat
dvec = array(0,N)
Amat = cbind(muvec, ones)
b0vec = c(Eport, 1)
if(noshort==1) {
idmat = matrix(0,N,N)
diag(idmat) = 1
Amat = cbind(Amat,idmat)
b0vec = c(b0vec, array(0,N))
}
wvec = solve.QP(Dmat, dvec, Amat, b0vec, meq=2)$solution
}
#initialize
N = 5
muvec = colMeans(returns[,1:N])
covmat = var(returns[,1:N])
#scroll through Eport values to derive efficient frontier
mincut = min(muvec)
maxcut = max(muvec)
Eportvec = seq(mincut,maxcut,length=300)
sigportvec = unlist(lapply(Eportvec, QPoptim, noshort=1, N, muvec, covmat))
#annualize stats
sigportvec = sigportvec * sqrt(12)
Eportvec = Eportvec * 12
#defining efficient frontier---------------------------------------------------
Emincut = Eportvec[which(sigportvec==min(sigportvec))]
idx = which(Eportvec>=Emincut)
#selecting the optimal portfolio
if (parameter=="mu") {
muportvec=c(d_mean)
w1 = meanvarweights(Eport=muportvec[1]/12, noshort=1, N, muvec, covmat)
sig1 = sqrt( t(w1) %*% covmat %*% w1 ) * sqrt(12)
} else {
sigportvec_1 = sigportvec[idx]
Eportvec_1 = Eportvec[idx]
a = which(abs(sigportvec_1-d_sd)==min(abs(sigportvec_1-d_sd)))
muportvec = c(Eportvec_1[a])
w1 = meanvarweights(Eport=(muportvec[1]/12), noshort=1, N, muvec, covmat)
sig1 = sqrt( t(w1) %*% covmat %*% w1 ) * sqrt(12)
}
#return data for output function-----------------------
temp = list(Eportvec = Eportvec, sigportvec = sigportvec,
w1=w1, muportvec=muportvec,sig1=sig1,idx=idx)
})
#####Plotting frontier output$prcPlot <- renderPlot({
#calls above function for prepped data------------------------------
temp = dataPrep()
sigportvec = temp$sigportvec
Eportvec = temp$Eportvec
muportvec = temp$muportvec
w1 = temp$w1
sig1 = temp$sig1
idx= temp$idx
#define variables from input boxes----------------------------------
startdt = ymd(input$start)
enddt = ymd(input$end)
parameter = input$parameter
d_mean = input$desired_annual_expected_return
d_sd = input$desired_annual_vol
#text heading of plot
startdt_txt = format(startdt, "%Y-%m-%d")
enddt_txt = format(enddt, "%Y-%m-%d")
str1 = "Efficient Frontier (based on data from)"
str2 = "to"
main_text_string = paste(str1,startdt_txt,str2,enddt_txt)
#plots all minimum variance portfolios------------------------------------------
plot(sigportvec, Eportvec,
xlim=c(0,max(sigportvec)), ylim=c(0,max(Eportvec)),
type="l", xlab="sigma", ylab="E(r)",
main=main_text_string, col = "black", lwd=1, lty="dashed")
#now just plot efficient frontier on top------------------------------------------
lines(x=sigportvec[idx], y=Eportvec[idx], type="l", col = "blue", lwd=2)
#point label text
sig1_d = format(round(sig1,2), nsmall = 2)
muportvec_d = format(round(muportvec,2), nsmall = 2)
sig1_t = as.character(sig1_d)
muportvec_t = as.character(muportvec_d)
sig1_txt = paste("sig=",sig1_t)
muportvec_txt = paste("mu=",muportvec_t)
label1 = paste(sig1_txt, ",", muportvec_txt)
#pick a portfolio along efficient frontier-----------------------
points(x=c(sig1), y=muportvec, col="blue", lwd=3, pch=1)
text(x=c(sig1), y=muportvec, labels = c(label1), pos=4)
})
#Making table heading output$titleTable <- renderTable({
#calls above function for prepped data------------------------------
temp = dataPrep()
sigportvec = temp$sigportvec
Eportvec = temp$Eportvec
muportvec = temp$muportvec
w1 = temp$w1
sig1 = temp$sig1
idx= temp$idx
#define variables from input boxes----------------------------------
startdt = ymd(input$start)
enddt = ymd(input$end)
parameter = input$parameter
d_mean = input$desired_annual_expected_return
d_sd = input$desired_annual_vol
#Prepare table heading format
sig1_d = format(round(sig1,2), nsmall = 2)
muportvec_d = format(round(muportvec,2), nsmall = 2)
if (parameter=="mu") {
string1 <- "The following portfolio achieves your desired annual mu ="
string2 <- "with vol ="
result = paste(string1, muportvec_d, string2, sig1_d)
} else {
string1 <- "The following portfolio achieves your desired annual vol ="
string2 <- "with mu ="
result = paste(string1, sig1_d, string2, muportvec_d)
}
#making table
colnamevec = c(" ")
numcol = length(colnamevec)
blank = array("", numcol)
blank = as.data.frame(t(blank))
colnames(blank) = colnamevec #prepared so that header of output is empty
row1 = blank
row1[1] = c(result)
ltemp = list(row1)
FINALOUT = rbindlist(ltemp) #prints this final table
}, align = 'c')
#####Plotting Summary Table output$prcTable <- renderTable({
#calls above function for prepped data------------------------------
temp = dataPrep()
sigportvec = temp$sigportvec
Eportvec = temp$Eportvec
muportvec = temp$muportvec
w1 = temp$w1
sig1 = temp$sig1
idx= temp$idx
#define variables from input boxes----------------------------------
startdt = ymd(input$start)
enddt = ymd(input$end)
parameter = input$parameter
d_mean = input$desired_annual_expected_return
d_sd = input$desired_annual_vol
#reducing decimal places in the weight values
w1_d = format(round(w1,2),nsmall=2)
#making table
colnamevec = c(" ", " ")
numcol = length(colnamevec)
blank = array("", numcol)
blank = as.data.frame(t(blank))
colnames(blank) = colnamevec #prepared so that header of output is empty
row1 = blank
row2 = blank
row3 = blank
row4 = blank
row5 = blank
row1[1:2] = c("MSFT",w1_d[1])
row2[1:2] = c("WMT",w1_d[2])
row3[1:2] = c("AAPL",w1_d[3])
row4[1:2] = c("IBM",w1_d[4])
row5[1:2] = c("KO",w1_d[5])
ltemp = list(row1, row2, row3, row4, row5)
FINALOUT = rbindlist(ltemp) #prints this final table
}, align = 'cc')
}
Run the application
shinyApp(ui = ui, server = server)
