It looks as if you're using the (discrete) inner product
$$\langle f,g\rangle=\sum_j^n f(x_j) g(x_j)$$
and I'll proceed with that assumption.
The canonical method for producing a basis of orthogonal polynomials with respect to some (discrete or continuous) inner product is the Stieltjes procedure. I discussed the continuous case here; the discrete case proceeds similarly.
In any event, to summarize: if you have a set of monic orthogonal polynomials $\phi_k(x)$ satisfying the recurrence
$$\phi_{k+1}(x)=(x+b_k)\phi_k(x)+c_k\phi_{k-1}(x)$$
with the associated orthogonality condition $\langle \phi_j,\phi_k\rangle=0,\quad j\neq k$ and the initial values $\phi_{-1}(x)=0$ and $\phi_0(x)=1$, the recurrence coefficients $b_k$ and $c_k$ are obtained through the formulae
$$\begin{align*}b_k&=-\frac{\langle x\phi_k ,\phi_k\rangle}{\langle \phi_k,\phi_k\rangle}\\c_k&=-\frac{\langle\phi_k,\phi_k\rangle}{\langle\phi_{k-1},\phi_{k-1}\rangle}\end{align*}$$
where $x\phi_k$ denotes the product of $x$ and $\phi_k(x)$.
For this particular case, as you say, $\phi_1(x)=x$ is easily determined; thus, to compute $\phi_2(x)$, we compute $b_1=0$ and $c_1=-\frac18$ to yield $\phi_2(x)=x^2-\frac18$.
We can then determine that the quadratic least-squares fit to the data is (approximately) given by $0.71+3.46\phi_1(x)+0.285714\phi_2(x)$.
If you want more details on data-fitting with orthogonal polynomials, see this article by George Forsythe. You might also want to look up the Gram polynomials, which are the standard orthogonal polynomials used for data fitting when the data has equispaced abscissas.
Best Answer
$AR(1)$ is one the most common models in longitudinal data, as unlike in time series where you have a lot realizations of the same process, in longitudinal data you usually have very few realization of multiple processes. Therefore, you don't have enough repeated measurements in order to perform valuable estimation of complex $ARMA(p,q)$ structures (large $p$ and $q$). In econometric analysis it us pretty common to find models like $$ y_{t+1}=\beta_0 + \beta_1y_{t} + \beta_2x+\epsilon_t, $$
that is, you have $AR(1)$ process as a part of the regression.