Is there a regression method for fitting general N-order polynomials of two or more variables? I would like to fit a general N-order polynomial of two or more variables. For example, consider (for reference) the general second-order polynomial in two variables $x$ and $y$:
$$\beta_{00}+\beta_{01}y^1+\beta_{02}y^2+\beta_{10}x+\beta_{11}xy+\beta_{12}xy^2+\beta_{20}x^2+\beta_{21}x^2y+\beta_{22}x^2y^2.$$
Another example could be the general N-order polynomial in two variables $x$ and $y$:
$$\sum_{n=0}^{N}\sum_{m=0}^{N}\beta_{nm}x^ny^m$$
I know that there is an extension of linear regression named Multiple Linear Regression, but is there a regression method for fitting general N-order polynomials of two or more variables?
 A: Consider the following two expressions:
$$z_i = \sum_{n=0}^{N}\sum_{m=0}^{N}\beta_{nm}x_i^ny_i^m+AE_i,\qquad i=1,...,T > N^2$$
where $AE_i$ stands for "approximation error",
and 
$$z_i = \sum_{n=0}^{N}\sum_{m=0}^{N}\beta_{nm}x_i^ny_i^m +u_i\;,\qquad E(u_i\mid X,Y) = 0\;,  i=1,...,T > N^2$$
where $u_i$ is a random variable with the stated properties, and $\{X,Y\}$ is to be understood as containing the series of $(x,y)$ but also of their powers etc.
If we write them down in full, in both cases we have a $T\times N^2$ system of equations that is linear in the $N^2$ unknown coefficients. 
What is their difference?
The first is NOT a stochastic specification, but a mathematical approximation specification: we want to express the values of $Z$ as a function of a $N$-order polynomial in $X,Y$. We do not make any assumptions as to whether any of $Z, X,Y$ are random or not, we just treat them as three $T\times 1$ column vectors containing numbers, vectors which we want to relate in the way the specification expresses.
The second is a stochastic specification, where even if we assume that the $X,Y$ are deterministic vectors, the existence of $U$ which is assumed random, makes $Z$ a random variable.
But in both cases, since a priori we have $T > N^2$, each system will be overidentified as a matter of basic mathematical results. This means that in either case we cannot find values for the unknown coefficients such that the vectors containing the terms $AE_i$ and $u_i$ respectively will be zero vectors.
Since we cannot achieve the ideal, it is natural to attempt to somehow optimize the way we will arrive at some values for the unknown coefficients.  
In the mathematical approximation specification, our motive is to achieve the "best" approximation, and that's it. Intuitively "best approximation" translates into "minimum approximation error". But trying to minimize the sum of $AE_i$'s would be misleading, since $AE_i$'s can be both positive and negative, and so, in a simple sum  they will tend to offset each other. The next natural step is to consider the sum of absolute values, which is indeed the most intuitive concept. But the absolute value is a nasty beast to manipulate mathematically, and this is how we arrived at the concept of minimizing the sum of squares -and the "least-squares" methodology was born (centuries ago). In the context of the mathematical approximation specification, this leads us to apply the "General Linear Least Squares"  which minimize $\sum_{i=1}^TAE^2_i$. Note again that the crucial matter in order to apply the least-squares methodology is whether the unknown coefficients (the subject of calculation/estimation) have a pure linear relationship with $Z$. We don't care whether the elements of $X,Y$ are squared, square-rooted or whatever.  
In the stochastic specification we postulate a deterministic relation between $Z$ on the one hand and $X,Y$ on the other, but one which is "disturbed" by unpredictable (hence "random") factors, embodied in $U$. We also assume (or hope), that "on average" these disturbances cancel out, hence the expectational assumption on $U$. So here, we assume that there exists a TRUE value of the coefficients beforehand describing the relationship, and we want to estimate this true value. Due to the contamination of the data by the disturbances, we cannot uncover the exact true values... the next best thing is for our estimates to be as close as possible to these true values... and one way to do that is to apply the least-squares methodology, which has strong intuitive power, in a "multiple linear regression" setting, as it is called in statistics.
To summarize: in the mathematical approximation specification, we care directly about the magnitude of the approximation error. In the stochastic specification, we care predominantly about the magnitude of the distance between the coefficient estimates and their true values, but we "discover" that this may be represented as isomorphic to the magnitude of "residual" which is the statistical term for the approximation error.
So we end up applying the same mathematical methodology on the same data, which means that the least-squares coefficient estimates in these two systems of equations will be exactly the same. 
But the goal and motive is different, and also, the subsequent interpretation of the results is different.
In the first case, the mathematical approximation is just that: we have described $Z$ as an approximate function of $X,Y$. Without further assumptions or arguments, we cannot use this results to say something about another data set containing values of $Z,X,Y$.
In the second case, the stochastic specification is endowed from the beginning with the assumption that there is a general relation between the $Z,X,Y$ variables, one that goes beyond the specific sets of numbers we have available. So the results we obtain can in prnciple be validly used to talk about this general relation that holds between these variables, and not just for the specific sample.
So if one starts from "I want to fit a $N$-order polynomial" one is closer to the spirit of the mathematical approximation method, because one does not argue as to why the $Z$ variable should somehow relate to the $X,Y$ set (including powers etc). We need a pre-estimation argument that involves either causality or at least association between $Z$ and $X,Y$ in order to be in the spirit of the stochastic specification.
A: The regression method for fitting general N-order polynomials of two or more variables is the General Linear Least Squares. Each monomial of the polynomial represents a basis function of the design matrix. For more details check General Linear Least Squares Technique
