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

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.

• Are you sure this is a constrained regression problem? As I read the question, you seek a relationship of the form $y_4$ (one Forex series) = $\beta_1 y_1 + \beta_2 y_2 + \beta_3 y_3$ (plus, I presume, a fourth term representing a prevailing safe rate of return). That's independent of the investment decision. If a customer wants to invest $c$ capital in $y_4$ using $y_1$, $y_2$, and $y_3$ as proxies, then they would just invest $c\beta_1$ in $y_1$, $c\beta_2$ in $y_2$, and $c\beta_3$ in $y_3$. That adds no special complication to the regression, does it? – whuber Jan 23 '12 at 17:03
• It does because if you model this you will find that B1 + B2 + B3 > 1 in many cases (or < 1 in others). That is because the currency one is trying to replicate with the descriptors will typically have a larger or smaller volatility than the others, and so the regression will give you smaller or larger weights in response. This requires the investor either not to be fully invested, or to leverage, which I do not want. As for safe rate of return no. All we are trying to do is replicate series1 using other variables. Being a finance guy and not a statistician perhaps I have misnamed my question. – Thomas Browne Jan 23 '12 at 18:49
• The reason for including a term for a safe rate of return is that sometimes it will have a nonzero coefficient. Presumably, safe instruments (overnight bank deposits) are available to everyone at low cost, so anyone ignoring this as a component of their investment basket could be choosing suboptimal combinations. Now, if the coefficients do not add to unity, so what? Just invest as much as you wish in the proportions estimated by the regression. – whuber Jan 23 '12 at 19:06
• right..... simple as that. Thanks. I feel a bit silly now haha. – Thomas Browne Jan 23 '12 at 19:11
• Not silly at all. Merely asking this question reflects a high level of thought. I was just checking my own understanding of your question to make sure you got an effective answer. Cheers. – whuber Jan 23 '12 at 19:14

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).

• Actually quick question. Shouldn't I be minimizing variance rather than sum? Isn't that what a regression does is minimize the variance of the square of the errors? – Thomas Browne Jan 23 '12 at 19:13
• This is clever, Elvis, but couldn't you accomplish the same thing simply by reparameterizing the regression? E.g., let $Y = \alpha_1 X_1 + \alpha_2 X_2 + (1-\alpha_1-\alpha_2)X_3 +\varepsilon$ That's equivalent to $Y-X_3 = \alpha_1(X_1-X_3) + \alpha_2(X_2-X_3)+\varepsilon$. The estimates and standard errors of the $\pi_i$ are straightforward to compute from the estimates and var-covar matrix of $\alpha_1$ and $\alpha_2$. – whuber Jan 23 '12 at 19:30
• @whuber Yes but with more noisy data, or with some of the $\pi_k$ close to $0$, you’d violate easily the constraint $\pi_k > 0$, which is the "hard" part of the problem. – Elvis Jan 23 '12 at 19:34
• A positive coefficient tells you to buy a foreign currency; a negative coefficient tells you to sell it. If you don't own that currency already, you need to borrow it in order to sell it ("selling short"). Because unrestricted borrowing can get people into trouble, there are constraints on the amount of borrowing and how it is paid for ("margin requirements" and "capital carrying costs" and "mark-to-market" procedures). Therefore, borrowing is possible but is often avoided except by major players in the markets or unless it confers large advantages. – whuber Jan 23 '12 at 19:50
• Many thanks to all for all the help. Actually just to make a comment on FX markets in general, they are more easy to short than equities or bonds because one does not have to borrow a stock before short selling. One simply flips the denominator and numerator currencies. So for example selling EURUSD and selling USDEUR are exactly equivalent trades in terms of the risk department, but they are of course exactly opposite positions. That's why FX is such a great playground for quant traders because they don't have to worry much about directional frictions which are much more important in equities – Thomas Browne Jan 24 '12 at 9:42

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))

• In the above case by Matifou, what's to prevent the third coefficient from being negative? For example, had the optimal coefficients for $\beta_1=0.75$ and $\beta_2=0.5$ we would get that $(1-\beta_1-\beta_2)=-0.25$ which implies here that our third coefficient is negative and therefore does not hold based on our desired regression. – A.S. Feb 21 '16 at 19:02
• Thanks @A.S. for pointing this out. Indeed, this solution works only for the equality constraints, not the inequality ones. I edited the text accordingly. – Matifou Feb 21 '16 at 20:30

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}$$

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.