Sum-of-squares decomposition and $R^2$ bounds

The two normal equations from the least-squares theorem above say $\sum_{i=1}^n e_i=0$ and $\sum_{i=1}^n x_ie_i=0$, where $e_i=y_i-\hat y_i$. Since $\hat y_i=\hat\beta_0+\hat\beta_1x_i$ is a linear combination of $1$ and $x_i$, both normal equations combine to give $\sum_{i=1}^n e_i\hat y_i=\hat\beta_0\sum_{i=1}^n e_i+\hat\beta_1\sum_{i=1}^n x_ie_i=0$. Combined with $\sum e_i=0$, this gives $\sum_{i=1}^n e_i(\hat y_i-\bar y)=\sum e_i\hat y_i-\bar y\sum e_i=0$.

Now expand $SST=\sum_{i=1}^n(y_i-\bar y)^2=\sum_{i=1}^n\big((y_i-\hat y_i)+(\hat y_i-\bar y)\big)^2=\sum e_i^2+2\sum e_i(\hat y_i-\bar y)+\sum(\hat y_i-\bar y)^2$. The cross term vanishes by the orthogonality just shown, leaving $SST=SSE+SSR$.

Since $SSE=\sum e_i^2\ge0$ and $SSR=\sum(\hat y_i-\bar y)^2\ge0$, dividing $SST=SSR+SSE$ by $SST>0$ gives $R^2=SSR/SST\in[0,1]$.

Finally, $\hat y_i-\bar y=\hat\beta_1(x_i-\bar x)$ (since $\hat y_i=\hat\beta_0+\hat\beta_1x_i$ and $\bar y=\hat\beta_0+\hat\beta_1\bar x$), so $SSR=\hat\beta_1^2S_{xx}$. Substituting $\hat\beta_1=S_{xy}/S_{xx}$ gives $SSR=S_{xy}^2/S_{xx}$, and since $SST=S_{yy}=\sum(y_i-\bar y)^2$, $R^2=\dfrac{S_{xy}^2}{S_{xx}S_{yy}}=\left(\dfrac{S_{xy}}{\sqrt{S_{xx}S_{yy}}}\right)^2=r^2$, the square of the sample correlation coefficient.