Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

This post explains how to use ROI and ROI.plugin.neos packages in R code, which provide an interface to NEOS. NEOS (Network-Enabled Optimization System) Server is a free internet-based service for solving numerical optimization problems. For understanding how to use ROI package, we implement R code for portfolio optimization problem.

# NEOS Optimization Solver in R code

NEOS can solve optimization problems with hundreds or thousands of control variables and constraints with so many time-tested algorithms. For more information, you can visit NEOS (https://neos-server.org/neos/).

To use NEOS service in R, ROI package and its plug in packages are needed. As we use ROI and ROI.plugin.neos package, install these packages in R studio.

### Point To Note

ROI provides too many functionalities to be covered in this post. It, also, has many plug-in packages as you can see the plug-in list from the above figure, Therefore this post deal with small part of it.

As our focus centers mainly on portfolio optimization or asset allocation. we consider a quadratic optimization problem.

For mean-variance portfolio optimization, we use Michaud dataset which are used on the previous post .

The mean-variance portfolio optimization problem has the following form.

\begin{align} & min_{\boldsymbol{w}} ⁡\frac{1}{2} \boldsymbol{w}^{‘}\Sigma \boldsymbol{w} \\ s.t. & \\ & \boldsymbol{w}^{‘} \boldsymbol{1} = 1 \\ & \boldsymbol{w}^{‘} \boldsymbol{\mu} = E(R) \\ & 0 \le \boldsymbol{w} \le 1 \end{align} $$\boldsymbol{w}$$ : weight vector
$$\Sigma$$ : covariance matrix
$$\mu$$ : mean return vector
$$E(R)$$ : target expected return

### ROI and ROI.plugin.neos

ROI package has the following format.

1) Objective Function

 123 # objective function to minimize    obj <– Q_objective(Q = 2*cov) Colored by Color Scripter cs

Objective function is set by using the function Q_objective() with arguments of Q or L or both. Q means the quadratic coefficient matrix and L means the linear coefficient vector.

2) Constraints

 123456 # 2 linear constraints    const <– L_constraint(         L   = rbind(rep(1, nvar), mu),         dir = c(“==”,“<="),         rhs = c(1, rset)) Colored by Color Scripter cs

Since portfolio optimization use the linear constraints, the function L_constraint() which has arguments of L, dir, and rhs is used. L means a left hand side coefficient matrix. dir means a vector of equalities or inequalities or both. rhs means a vector of right hand side values. If quadratic constraints are used, function Q_constraint() is used.

3) Upper and Lower Bounds

 1234 # lower and upper bounds    lb_ub <– V_bound(lb = rep(0, nvar),                     ub = rep(1, nvar)) Colored by Color Scripter cs

Bounds or ranges for control variables are set by using the function V_bound() which has two arguments, lb and ub. This usage is straightforward. lb and ub mean vectors of lower and upper bounds for control variables respectively.

4) Optimization Problem

 12345 # create optimization problem    op <– OP( objective   = obj,               constraints = const,               bounds      = lb_ub) Colored by Color Scripter cs

Final optimization problem is set by using the function OP(). This function take the above three arguments as input arguments.

5) Solving Problems using NEOS Solver

 123456 res_scip  <– ROI_solve(op, solver = “neos”, method = “scip”,                            email=“your email address”)    res_mosek <– ROI_solve(op, solver = “neos”, method = “mosek”,                           email=“your email address”)       res_cplex <– ROI_solve(op, solver = “neos”, method = “cplex “,                           email=“your email address”) cs

Using function ROI_solve() with solver = “neos”, we can use NEOS solver service. Of course, arguments such as solver and method have several options. In particular, method can take several optimizers. In this case, we use three alternative optimizers: scip, mosek, cplex. These are time-tested optimizers for the large-scale optimization problems. Here, solver = “neos” requires the package ROI.plugin.neos.

From recently, user email address is required so that you sign up and log in the NEOS website and your email should be verified. Every time you run ROI_solve() with solver = “neos”, You receive email which contains the full results from NEOS

### R code using NEOS solver through ROI

The following R code is implemented based on the above format. Of course, traditional optimization using function solveQPXT() is the same as the previous post regarding Michaud portfolio optimization.

 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117 #=================================================================## Finance and Insurance Engineering using R # by Sang-Heon Lee & Hosam Ki## https://kiandlee.blogspot.com#—————————————————————–## Portfolio Optimization using ROI NEOS#=================================================================# graphics.off()  # clear all graphsrm(list = ls()) # remove all files from your workspace library(quadprogXT)      # solveQPXTlibrary(ROI)             # ROI_solvelibrary(ROI.plugin.neos) # NEOS #—————————————————–# Michuad dataset   #—————————————————–     setwd(“D:/SHLEE/blog/rneos/michuad”)     # dataset in Michaud and Michaud (2007) appendix    df.data <– read.csv(‘mu_sd_corr_michaud.csv’)    head(df.data)        # mean, std, correlation    mu   <– as.vector(df.data[,1])    sd   <– as.vector(df.data[,2])    corr <– as.matrix(df.data[,–c(1,2)])        # number of asset    nvar <– length(mu)    var.name <– colnames(df.data[,–c(1,2)])        # convert correlation to covariance matrix    cov <– diag(sd)%*%corr%*%diag(sd) #—————————————————–# Portfolio optimization using solveQPXT#—————————————————–        n.er <– 100  # number of EF points    rset <– seq(min(mu),max(mu),length=n.er+2)    rset <– rset[2:n.er+1]        # given returns and unknown std and weight    port1.ret <– rset    port1.std <– rset*0    port1.wgt <– matrix(0,n.er,nvar)        # ith portfolio problem setting    i = 10;    Dmat <– 2*cov    dvec <– rep(0,nvar) #c(0,0)    Amat <– t(rbind(t(rep(1,nvar)),t(mu),diag(nvar)))    bvec <– c(1,rset[i],rep(0,nvar))            # mean-variance optimization    m<–solveQPXT(Dmat,dvec,Amat,bvec,meq=2,factorized=FALSE)            # output    port1.std[i]  <– sqrt(m$value) port1.wgt[i,] <– t(m$solution)    #—————————————————–# Portfolio optimization using ROI.neos.plugin#—————————————————–          #——————–    # min w’*cov*w      # s.t.     #    sum(w) = 1    #    w*r   <= mu    #    0 <= w <= 1    #——————–        # objective function to minimize    obj <– Q_objective(Q = 2*cov)           # 2 linear constraints    const <– L_constraint(         L   = rbind(rep(1, nvar), mu),         dir = c(“==”,“<="),         rhs = c(1, rset))            # lower and upper bounds    lb_ub <– V_bound(lb = rep(0, nvar),                     ub = rep(1, nvar))        # create optimization problem    op <– OP( objective   = obj,               constraints = const,               bounds      = lb_ub)            res_scip  <– ROI_solve(op, solver = “neos”, method = “scip”,                            email=“your email address”)    res_mosek <– ROI_solve(op, solver = “neos”, method = “mosek”,                           email=“your email address”)       res_cplex <– ROI_solve(op, solver = “neos”, method = “cplex “,                           email=“your email address”)   #—————————————————–# Comparisons#—————————————————–        # function value    cat(paste(“\nobj val : solveQPXT   =”, m$value, “\n”, ” NEOS(scip) =”, res_scip$objval, “\n”,              ”         NEOS(mosek) =”, res_mosek$objval, “\n”, ” NEOS(cplex) =”, res_cplex$objval, “\n”))                    # weight comparisons : solveQPXT versus NEOS(mosek)    output_mosek <– cbind(res_mosek$solution, port1.wgt[i,], res_mosek$solution – port1.wgt[i,])    colnames(output_mosek) <– c(“NEOS(mosek)”, “solveQPXT”, “diff”)    print(output_mosek)    Colored by Color Scripter cs

Seeing the results below, we can find that although MOSEK optimizer attains the most minimized function value, the differences between optimizers are negligible. Therefore the optimal weight differences between solveQPXT and MOSEK is also so small and insignificant.

Too small weight differences can also be found the following bar plot.

### Final Thoughts

From this post, we can learn how to use NEOS optimizer using ROI and ROI.plugin.neos packages.

It is worth noting that NEOS is not for this small problem but for so big problem which has so many variables and constraints. Therefore, when you have to solve this big problem, it is suitable to use NEOS solver using R and its packages.

Before anything else, this NEOS solver can be used to the multi-stage stochastic linear programming because it can solve the problems written by GAMS or AMPL format using MOSEK or CPLEX. $$\blacksquare$$