Setting up arbitrage strategy in R

Question

I am trying to construct an arbitrage portfolio $\textbf{x}$ such that $S^T\textbf{x} = 0$ and $A\textbf{x} \geq \textbf{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. I am not sure how to code it up myself either. Any help or hints on how to proceed would be much appreciated. Thanks!

Answer 1

A test for arbitrage opportunities with an LP is to minimize the cost of setting up the portfolio, subject to the restriction that the portfolio loses money in no state of the world. (Note that in your formulation you are missing the actual objective; you only list constraints.) If you find a portfolio that has a negative cost (i.e. you get paid for holding it), but you never lose money, you have found an arbitrage portfolio. Or if find a portfolio with zero cost, but no possibility of losses and at least one positive payoff, you have found an arbitrage opportunity. If you find one arbitrage portfolio, then without constraints you'll typically find infinitely many. That should be intuitive: if you have a portfolio of zero cost, but only non-negative payoffs, you can multiply all weights by some constant and still have an arbitrage portfolio. Also, without constraints, an arbitrage portfolio with negative cost will be unbounded.

Doing this in R:

A <- matrix(c(2, 1, 0, 3, 1,
              1, 1, 1, 2, 1,
              0, 1, 2, 1, 0), byrow = TRUE, nrow = 3)
S <- c(1, 1, 1, 2, 1/3)


library("Rglpk")
bounds <- list(lower = list(ind = 1:5, val = rep(-Inf, 5)))
lp.sol <- Rglpk_solve_LP(S,
                         mat = A,
                         dir = rep(">=", 3),
                         rhs = c(0, 0, 0),
                         bounds = bounds,
                         control = list(canonicalize_status = FALSE,
                                        verbose = TRUE))
## [....]
## LP HAS UNBOUNDED PRIMAL SOLUTION

That's not overly helpful because it only tells you that there is an arbitrage opportunity. So we add constraints: a negative position may not exceed -1.

bounds <- list(lower = list(ind = 1:5, val = rep(-1, 5)))
lp.sol <- Rglpk_solve_LP(S,
                         mat = A,
                         dir = rep(">=", 3),
                         rhs = c(0, 0, 0),
                         bounds = bounds)
sum(lp.sol$solution*S)
## [1] -1
A %*% lp.sol$solution
##      [,1]
## [1,]    0
## [2,]    3
## [3,]    0

Now you have a negative-cost portfolio (i.e. you receive 1 for setting up the portfolio). To make it zero cost, you invest those proceeds into one asset:

x <- lp.sol$solution
x[2] <- x[2] + 1/S[2]
sum(x*S)
## [1] 0
A %*% x
##      [,1]
## [1,]    1
## [2,]    4
## [3,]    1

Now you have a zero-cost portfolio with strictly positive payoffs.

Alternatively, you could use another numerical solver to directly solve the optimization model. Here is an example. (Disclosure: I am the maintainer of packages NMOF and neighbours.) It's more convenient to work with returns:

R <- t(t(A)/S) - 1
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0   -1  0.5    2
## [2,]    0    0    0  0.0    2
## [3,]   -1    0    1 -0.5   -1


library("NMOF")         ## https://github.com/enricoschumann/NMOF
library("neighbours")   ## https://github.com/enricoschumann/neighbours

Now we directly maximize the average payoff, say. (The implementation I use minimizes, so I multiply by -1.)

max_payoff <- function(x, R, S)
    -sum(R %*% x) +                ## => maximize average payoff
    -10*sum(pmin(R %*% x, 0))      ## => penalty for negative state returns

nb <- neighbourfun(-1, 5, length = 5, stepsize = 5/100)

ta.sol <- LSopt(max_payoff,
                list(neighbour = nb,
                     x0 = rep(0, length(S)),
                     nI = 5000),
                R = R, S = S)
round(ta.sol$xbest, 3)  ## the portfolio
## [1] -1.00 -1.00  0.75 -1.00  2.25
round(R %*% ta.sol$xbest, 1)  ## the state returns
##      [,1]
## [1,]  2.2
## [2,]  4.5
## [3,]  0.0

The portfolio in shares:

x <- round(ta.sol$xbest/S, 3)
sum(x*S)
## [1] 0
A %*% x
##      [,1]
## [1,] 2.25
## [2,] 4.50
## [3,] 0.00