Optimal polynomial order in equality-constrained linear regression In fitting a polynomial $p_N = \sum\limits_i^N a_i x^i$ to my dataset $(x,y)$ using a least squares criteria, I'm interested in determining the optimal polynomial order $N$. The caveat is that, owing to problem specific restrictions, the polynomial must satisfy the following boundary conditions: $p_N(0) = 0$, and $\frac{d}{dx}p_N(L)=0$. ($L$ is known).
I gather that if the problem was unconstrained, the order could be determined using 2 successive ANOVA F-tests: the first to determine whether or not order $N$ adequately explains the variance in the observed $Y$, and if it passes, a second test to determine whether the added explained variance by order $N$ over order $N-1$ is statistically significant.
Is there a similar test that I could apply that accounts for the presence of boundary conditions? 
Thank you in advance for your time!
 A: As a warning - fitting polynomials like this is fraught with danger in several respects.  Even for smallish $N$ you risk having an ill-conditioned problem, especially if the range of $x$ values is large, with substantial, perhaps total, loss of significant digits of accuracy.  Additionally, polynomials often are poor functions for approximating functions that are not themselves polynomials.  You often wind up with way too much "wiggliness" in regions where there aren't many x-values.  Splines tend to do much better in both regards (as in R package mgcv), and also get around the problem of selecting polynomial degree.  Unfortunately, automated spline estimation packages won't help you with your derivative constraint, at least not easily.
Your testing procedure also relies on an $F-$test to determine whether a polynomial of degree $N$ adequately explains the variance in $y$, but that won't work.  To see this, consider $y = x + e$, a first degree polynomial, with $x$ distributed, say, $U(0,1)$.  The variance of $y$ depends very heavily on the variance of $e$, in fact, making standard assumptions, $\sigma^2_y = \sigma^2_e + 1/12$.  The amount of the variance of $y$ that can be explained by the polynomial in $x$ can be arbitrarily close to either 0% or 100%, depending on $\sigma^2_e$, but whatever that amount is, the model $y = x + e$ is correct.  Clearly no test, other than human judgment, can determine whether the amount of the variance of $y$ that can be explained by the polynomial in $x$ is adequate or not.  The issue of the definition of adequate in a mathematically formal way also arises; adequate is determined by some external, human, objective, and is often not readily quantifiable.
On to your question - the constrained polynomial equation can be transformed into an unconstrained equation with two fewer parameters, then standard techniques can be applied to the resultant equation.  Consider the boundary constraint $P_N(0)=0$.  Your expression has a constant term; just eliminate the constant term from the polynomial, and, when $x=0$, the polynomial will equal 0 as well.  The derivative constraint can be incorporated directly into the expression by writing out the derivative and observing that the parameter $a_1$ can be expressed as a polynomial in $L$:  $a_1 = - \sum_{i=2}^Nia_iL^{i-1}$.  This relationship has to hold for the derivative constraint to be satisfied, and vice versa.  Then you can substitute this expression for $a_1$, rearrange terms, and estimate the resulting unconstrained linear equation instead.  For example, for $N=3$:  
$a_1 = -2a_2L -3a_3L^2$
$p_N = (-2a_2L -3a_3L^2)x + a_2x^2 + a_3x^3 = (-2Lx+x^2)a_2 + (-3L^2x + x^3)a_3$
By creating new right hand side variables, e.g., $z_1 = -2Lx+x^2$, you recover the linear regression form $p_N = a_2z_1 + a_3z_2$ etc., and you can run that regression instead, testing for improvement as you increment $N$.
