How do I fit a constrained regression in R so that coefficients total = 1?

Question

I see a similar constrained regression here:

Constrained linear regression through a specified point

but my requirement is slightly different. I need the coefficients to add up to 1. Specifically I am regressing the returns of 1 foreign exchange series against 3 other foreign exchange series, so that investors may replace their exposure to that series with a combination of exposure to the other 3, but their cash outlay must not change, and preferably (but this is not mandatory), the coefficients should be positive.

I have tried to search for constrained regression in R and Google but with little luck.

Answer 1

If I understand correctly, your model is $$ Y = \pi_1 X_1 + \pi_2 X_2 + \pi_3 X_3 + \varepsilon, $$ with $\sum_k \pi_k = 1$ and $\pi_k\ge0$. You need to minimize $$\sum_i \left(Y_i - (\pi_1 X_{i1} + \pi_2 X_{i2} + \pi_3 X_{i3}) \right)^2 $$ subject to these constraints. This kind of problem is known as quadratic programming.

Here a few line of R codes giving a possible solution ($X_1, X_2, X_3$ are the columns of X, the true values of the $\pi_k$ are 0.2, 0.3 and 0.5).

library("quadprog");
X <- matrix(runif(300), ncol=3)
Y <- X %*% c(0.2,0.3,0.5) + rnorm(100, sd=0.2)
Rinv <- solve(chol(t(X) %*% X));
C <- cbind(rep(1,3), diag(3))
b <- c(1,rep(0,3))
d <- t(Y) %*% X  
solve.QP(Dmat = Rinv, factorized = TRUE, dvec = d, Amat = C, bvec = b, meq = 1)
$solution
[1] 0.2049587 0.3098867 0.4851546

$value
[1] -16.0402

$unconstrained.solution
[1] 0.2295507 0.3217405 0.5002459

$iterations
[1] 2 0

$Lagrangian
[1] 1.454517 0.000000 0.000000 0.000000

$iact
[1] 1

I don’t know any results on the asymptotic distribution of the estimators, etc. If someone has pointers, I’ll be curious to get some (if you wish I can open a new question on this).

Answer 2

As mentioned by whuber, if you are interested only in the equality constraints, you can also just use the standard lm() function by rewriting your model:

\begin{eqnarray} Y&=&\alpha+\beta_1 X_1+\beta_2 X_2+\beta_3 X_3+\epsilon\\ &=& \alpha+\beta_1 X_1+\beta_2 X_2+(1-\beta_1-\beta_2) X_3+\epsilon\\ &=& \alpha + \beta_1( X_1-X_3) +\beta_2 (X_2-X_3)+ X_3+\epsilon \end{eqnarray}

But this does not guarantee that your inequality constraints are satisfied! In this case, it is however, so you get exactly the same result as using the quadratic programming example above (putting the X3 on the left):

X <- matrix(runif(300), ncol=3)
Y <- X %*% c(0.2,0.3,0.5) + rnorm(100, sd=0.2)
X1 <- X[,1]; X2 <-X[,2]; X3 <- X[,3]
lm(Y-X3~-1+I(X1-X3)+I(X2-X3))

Answer 3

Old question but since I'm facing the same problem I thought to post my 2p...

Use quadratic programming as suggested by @Elvis but using sqlincon from the pracma package. I think the advantage over quadrpog::solve.QP is a simpler user interface to specify the constraints. (In fact, lsqlincon is a wrapper around solve.QP).

Example:

library(pracma)

set.seed(1234)

# Test data
X <- matrix(runif(300), ncol=3)
Y <- X %*% c(0.2, 0.3, 0.5) + rnorm(100, sd=0.2)

# Equality constraint: We want the sum of the coefficients to be 1.
# I.e. Aeq x == beq  
Aeq <- matrix(rep(1, ncol(X)), nrow= 1)
beq <- c(1)

# Lower and upper bounds of the parameters, i.e [0, 1]
lb <- rep(0, ncol(X))
ub <- rep(1, ncol(X))

# And solve:
lsqlincon(X, Y, Aeq= Aeq, beq= beq, lb= lb, ub= ub)

[1] 0.1583139 0.3304708 0.5112153

Same results as Elvis's:

library(quadprog)
Rinv <- solve(chol(t(X) %*% X));
C <- cbind(rep(1,3), diag(3))
b <- c(1,rep(0,3))
d <- t(Y) %*% X  
solve.QP(Dmat = Rinv, factorized = TRUE, dvec = d, Amat = C, bvec = b, meq = 1)$solution

EDIT To try to address gung's comment here's some explanation. sqlincon emulates matlab's lsqlin which has a nice help page. Here's the relevant bits with some (minor) edits of mine:

X Multiplier matrix, specified as a matrix of doubles. C represents the multiplier of the solution x in the expression C*x - Y. C is M-by-N, where M is the number of equations, and N is the number of elements of x.

Y Constant vector, specified as a vector of doubles. Y represents the additive constant term in the expression C*x - Y. Y is M-by-1, where M is the number of equations.

Aeq: Linear equality constraint matrix, specified as a matrix of doubles. Aeq represents the linear coefficients in the constraints Aeq*x = beq. Aeq has size Meq-by-N, where Meq is the number of constraints and N is the number of elements of x

beq Linear equality constraint vector, specified as a vector of doubles. beq represents the constant vector in the constraints Aeq*x = beq. beq has length Meq, where Aeq is Meq-by-N.

lb Lower bounds, specified as a vector of doubles. lb represents the lower bounds elementwise in lb ≤ x ≤ ub.

ub Upper bounds, specified as a vector of doubles. ub represents the upper bounds elementwise in lb ≤ x ≤ ub.

Answer 4

As I understand your model, you're seeking to find $$ \bar{\bar{x}} \cdot \bar{b} = \bar{y} $$ such that $$ \sum \left [ \begin{matrix} \bar{b} \end{matrix} \right ] =1 $$

I've found the easiest way to treat these sorts of problems is to use matrices' associative properties to treat $\bar{b}$ as a function of other variables.

E.g. $\bar{b}$ is a function of $\bar{c}$ via the transform block $\bar{\bar{T_c}}$. In your case, $r$ below is $1$. $$ \bar{b} = \left [ \begin{matrix} k_0 \\ k_1 \\ k_2 \end{matrix} \right ] = \bar{\bar{T_c}} \cdot \bar{c} = \left [ \begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -1 & -1 & 1 \end{matrix} \right ] \cdot \left[ \begin{matrix} k_0 \\ k_1 \\ r \end{matrix} \right ] $$ Here we can separate our $k$nowns and $u$nknowns. $$ \bar{c} = \left[ \begin{matrix} k_0 \\ k_1 \\ r \end{matrix} \right ] = \bar{\bar{S_u}} \cdot \bar{c_u} + \bar{\bar{S_k}} \cdot \bar{c_k} = \left[ \begin{matrix} 1 & 0 \\ 0 & 1 \\ 0 & 0 \end{matrix} \right] \cdot \left [ \begin{matrix} k_0 \\ k_1 \end{matrix} \right ] + \left [ \begin{matrix} 0 \\ 0 \\ 1 \end{matrix} \right ] \cdot r $$ While I could combine the different transform/separation blocks, that gets cumbersome with more intricate models. These blocks allow knowns and unknowns to be separated. $$ \bar{\bar{x}} \cdot \bar{\bar{T_c}} \cdot \left ( \bar{\bar{S_u}} \cdot \bar{c_u} + \bar{\bar{S_k}} \cdot \bar{c_k} \right ) = \bar{y} \\ \bar{\bar{v}} = \bar{\bar{x}} \cdot \bar{\bar{T_c}} \cdot \bar{\bar{S_u}} \\ \bar{w} = \bar{y} - \bar{\bar{x}} \cdot \bar{\bar{T_c}} \cdot \bar{\bar{S_k}} \cdot \bar{c_k} $$ Finally the problem is in a familiar form. $$ \bar{\bar{v}} \cdot \bar{c_u} = \bar{w} $$

Answer 5

Using matrix algebra it is possible write following formula if you want to relax non-negative coefficients constraint,

$\beta=(X^{T}X)^{-1}X^{T}y+1\left[\frac{1_{scalar}-1^{T}(X^{T}X)^{-1}X^{T}y}{1^{T}(X^{T}X)^{-1}1}\right](X^{T}X)^{-1}.$

This might be helpful in case of a need for quick, simple and exact solution.

Answer 6

Plenty of different solutions already, but here's another that may be of use.

We can define a regression of the form $$Y \sim \text{Normal}(\mu,\sigma) \\ \mu = X \cdot \text{softmax}(\beta)$$

Applying a softmax function over the $\beta$ parameters constrains these to sum to one. The model is only a slight modification of a standard linear regression, which we can fit using MLE with optim.

set.seed(1234)

# Simulate
X <- matrix(runif(300), ncol=3)
y <- X %*% c(0.2,0.3,0.5) + rnorm(100, sd=0.2)

# Simple softmax function
softmax <- function(x) exp(x) / sum(exp(x))

# Normal linear regression log-likelihood, with softmax on betas
ll <- function(par, X, y){
    beta <- par[-1]
    sigma <- sqrt(par[1])
    mu <- X %*% softmax(beta)
    -sum(dnorm(y, mean = mu, sd = sigma, log = TRUE))
}

# Fit it
fit <- optim(
    par = rep(1, 4), 
    fn = ll,
    X = X,
    y = y
)

The $\beta$ parameters returned are unconstrained, but we can get the constrained estimates by transforming them using the softmax function (ignoring the unconstrained $\sigma$ term)

round(softmax(fit$par[-1]), 5)
[1] 0.15833 0.33047 0.51120

These results are quite consistent with Elvis and Dariober's quadratic programming methods above.

This was inspired by a blog post from Peter Ellis who discussed several of the methods in this thread, and introduced a clever Bayesian modelling solution using Dirichlet priors.