Search code examples
rfinance

Setting up an arbitrage strategy in R by solving system of inequalities subject to an equality constraint


I am trying to construct an arbitrage portfolio x such that Sx = 0 and Ax>=0, where A is the payoff matrix at t=1 and S is the price at t=0. I was not able to do it manually, so I tried using functions contained in the limSolve and lpSolve packages in R with no success, as I keep getting the zero-vector (I need nontrivial solutions). I am not sure how to code it up myself either. Any help or hints on how to proceed would be much appreciated. Thanks!

    A = data.frame(
              cbind(
                 c(2,1,0),
                 c(1,1,1),
             c(0,1,2),
             c(3,2,1),
             c(1,1,0)
              )
             ) %>% as.matrix()
f.con = A
    
    S = data.frame(
              cbind(
             c(1,1,1,2,1/3)
              )
             ) %>% as.matrix()

f.obj = c(t(S))

# Set unequality signs
f.dir <- c(">",
           ">",
           ">")

# Set right hand side coefficients
f.rhs <- c(0,
           0,
           0)

# Final value 
lp("min", f.obj, f.con, f.dir, f.rhs)$solution

# Variables final values
lp("max", f.obj, f.con, f.dir, f.rhs)$solution

Solution

  • We first explain why the code in the question got the results that it did.

    Note that there is an implicit constraint of x >= 0.

    Regarding the subject that refers to a equality constraint there are no equality constraints in the code shown in the question.

    Regarding the minimum, in the lp arguments > means >= so clearly x=0 is feasible. Given that all components of the objective vector are positive it leads to a minimum of 0. From ?lp

    const.dir: Vector of character strings giving the direction of the constraint: each value should be one of "<," "<=," "=," "==," ">," or ">=". (In each pair the two values are identical.)

    Regarding the maximum, there are no upper bound constraints on the solution and the objective vector's components are all positive so there is no finite solution. Whether the linear program was successful or not should always be shown prior to displaying the solution and if it was unsuccessful then the solution should not be displayed since it has no meaning.

    Regarding the code, cbind already produces a matrix so it is pointless to convert it to a data frame and then back to a matrix. Also the objective can be expressed as a plain vector and the rhs of the constraints can be written as a scalar which will be recycled to the appropriate length. We can write the code in the question equivalently as:

    library(lpSolve)
    
    A <- cbind(2:0, 1, 0:2, 3:1, c(1, 1, 0))
    S <- c(1, 1, 1, 2, 1/3)
    
    res1 <- lp("min", S, A, ">=", 0)   
    res1
    ## Success: the objective function is 0 
    res1$solution
    ## [1] 0 0 0 0 0
    
    res2 <- lp("max", S, A, ">=", 0)
    res2
    ## Error: status 3 
    

    CVXR

    It is easier to formulate this using CVXR as shown below.

    Find a vector x which satisfies Ax >= 0 and Sx == 0. (A and S are from above.) We add constraints -1 <= x <= 1 to keep the solution within bounds and use an arbitrary objective sum(x) since we are only looking for any feasible solution.

    library(CVXR)
    
    x <- Variable(5)
    objective <- Maximize(sum(x))  # arbitrary objective
    constraints <- list(A %*% x >= 0, sum(S*x) == 0, x >= -1, x <= 1)
    problem <- Problem(objective, constraints)
    soln <- solve(problem)
    

    giving:

    > soln$status
    [1] "optimal"
    
    > soln$value
    [1] 1.66666
    
    > soln$getValue(x)
               [,1]
    [1,]  0.8788689
    [2,]  0.4790521
    [3,]  0.3087133
    [4,] -0.9999857
    [5,]  1.0000113
    

    lpSolve again

    To do this using lpSolve make the change of variables

    x = 2*y-1
    

    or equivalently

    y = (x+1)/2
    

    which converts

    Ax >= 0
    Sx == 0
    -1 <= x <= 1
    

    to

    2Ay >= A1
    2Sy >= S'1
    0 <= y <= 1
    

    so we write:

    f.con <- rbind(2*A, 2*S, diag(5))
    f.dir <- c(rep(">=",3),  "=", rep("<=", 5))
    f.rhs <- c(rowSums(A), sum(S), rep(1, 5))
    
    res3 <- lp("max", rep(1, 5), f.con, f.dir, f.rhs)
    res3
    ## Success: the objective function is 3.333333 
    2*res3$solution-1
    ## [1]  1.000000e+00 -4.966028e-13  6.666667e-01 -1.000000e+00  1.000000e+00