1
$\begingroup$

I have a (slightly simplified) model of the following form:

$Y=c_1X_1 + c_2X_2 + \epsilon$ subject to the constraint $0\leq c_1\leq c_2$.

The distribution of $\epsilon$ is actually not important to me - fitting the curve by least squares is all I care about (for a volatility surface parameterization).

What is the best approach to solving this with deterministic results (i.e. not using iterative techniques dependent on starting points etc.).

I am working in C++ and I have Armadillo for my linear algebra library.

I am comfortable with python's libraries as well and R if necessary, but ultimate implementation will be C++.

$\endgroup$
7
  • 1
    $\begingroup$ Welcome to our site. As described in the answer at stats.stackexchange.com/questions/29279, write your model as $$Y = c_1(X_1 + X_2) + (c_2-c_1)X_2 + \epsilon$$ subject to the constraint $(c_1,c_2-c_1)\ge(0,0)$ and search the solutions for "non-negative least squares" problems given at stats.stackexchange.com/search?q=non-negative+least+squares. $\endgroup$
    – whuber
    Aug 28 '18 at 18:46
  • $\begingroup$ Thx - I knew it had to be something easy like this. I'll be back in the swing of things soon - been a while since doing this kind of stuff. I have to solve this problem a few thousand time a second, so we shall see how that goes. len(Y) usually between 10-100. If speed is a problem, I can do usual OLS and it would be ok if sometimes the constraint fails as long enough of them do not - usually that is the case. $\endgroup$ Aug 28 '18 at 18:59
  • $\begingroup$ When the OLS solution doesn't work, you know at least one constraint holds. It isn't too inefficient just to solve the resulting three possible models $Y=c_1(X_1+X_2)+\epsilon,$ $Y=c_2 X_2 + \epsilon,$ and $Y=\epsilon$ and select the one with least sum of squared residuals. You could probably do up to a million of these a second (per core) by coding the procedure in a compiled language and perhaps a billion a second using a couple of high-end coprocessors, assuming you have the bandwidth to get all the data in quickly enough! $\endgroup$
    – whuber
    Aug 28 '18 at 19:16
  • $\begingroup$ Good idea - when my data is good, the $R^2$ does need to be better than 0.99, so what I do is test a model such as $Y=c_1\rho X_1 + c_1 X_2$ where $|\rho|\leq 1$. I test for 20 or 40 different $\rho$ from -1 to 1. (Note that this like solving the original problem where $0\leq |c_1| \leq c_2$ which is what I am really doing. $\endgroup$ Aug 28 '18 at 19:27
  • 1
    $\begingroup$ That looks like a hard way to do it. Why not just use OLS to find $c_1$ and $c_2$ directly? $\endgroup$
    – whuber
    Aug 28 '18 at 19:28
1
$\begingroup$

Reformulate your model to $$ Y = c_1(X_1+X_2) + (c_2-c_1) X_2 + \epsilon $$ (and probably an intercept which you have left out). Then the constraints is $0\le c_1, 0\le (c_2-c_1)$, and you can use an implementation of nonnegative least squares. In R there is a package colf, see an example in Looking for function to fit sigmoid-like curve. See also Linear Regression with individual constraints in R.

(This approach will be much simpler than what you propose in comments.)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.