Imposing non-negativity constraint on a linear regression function

Question

Suppose I am interested in estimating the linear regression model

$$ Y_i = g(X_i)^T\beta + \epsilon_i $$

where $Y_i$ is a scalar outcome of interest, $X_i$ is a scalar covariate with support on the unit interval, $g(\cdot)$ is a $K$-dimensional vector of known functions that are not perfectly colinear, $\beta$ is a $K$-dimensional vector of parameters to be estimated, and $E(g(X_i) \epsilon_i) = 0$. Suppose I know that $\Pr(Y \geq 0) = 1$, so I'd like to impose the condition

$$ \beta^T g(x) \geq 0 \qquad \text{for all $x \in [0,1]$} $$

when estimating $\beta$. How would I go about doing this?

Answer 1

You may wish to solve the constrained Maximum likelihood problem: $$ \min_{\beta} \sum_{i=1}^n{(y_i-g(x_i)^T\beta)^2} \;\;\;s.t.\;\; \beta^Tg(x)\geq 0. $$ In the case that the dimensionality $K$ is fixed, the above estimator should be asymptotically efficient. How such estimators can be computed, of course, depend on the function $g$ and most importantly whether $\{\beta:g(x)^T\beta\geq 0\}$ is convex.

Answer 2

Assuming $g$ is polynomial, one way to attack the functional constraint on non-negativity is to employ a sum-of-squares approach. A sum-of-squares approach (conservatively) replaces a non-negativity constraint $g(x)\geq 0$ with a condition that the polynomial is a sum of squares. When non-negativity only is required on a region $h(x) \geq 0$, it can be generalized by the positivstellensatz to the search of a non-negative polynomial $s(x)$ such that $g(x)\geq s(x)h(x)$, once again replacing non-negativity of $s(x)$ with a sum-of-squares condition.

Sum-of-squares decompositions and methods based on this is an active area with publically available software. The code below implements this in the MATLAB toolbox YALMIP (disclaimer, developed by me)

% Define some data for a noisy quartic
xdata = 0:0.05:1;
Y = 5*(xdata-.5).^4 + .02*randn(1,length(xdata));
clf;plot(xdata,Y);grid

% Create monomial data
X = [xdata.^0;xdata;xdata.^2;xdata.^3;xdata.^4];

% Define decision variables and create residuals
beta = sdpvar(5,1);
e = Y - beta'*X;

% Least squares solution
optimize([],e*e');
hold on;plot(xdata,value(beta)'*X);

% Parameterize a quartic polynomial
sdpvar x
g = beta'*[1;x;x^2;x^3;x^4];

% g should be non-negative when  (x-0.5)^.2 <= 0.25
% define multiplier for positivstellensatz
[s,coeffs] = polynomial(x,2);

solvesos([sos(g - s*(.25-(x-.5)^2)),sos(s)],e*e',[],[beta;coeffs])

hold on;plot(xdata,value(beta)'*X);

You can read more here

https://yalmip.github.io/tutorial/sumofsquaresprogramming/