[proofplan]
We center the two random variables so that the variance of $Y_\beta$ becomes the second moment of a linear combination of centered variables. Expanding the square and using the definitions of variance and covariance gives a real quadratic polynomial in $\beta$. Since the coefficient of $\beta^2$ is $\operatorname{Var}(C)>0$, we minimize the quadratic by completing the square. Substituting the minimizing coefficient gives the stated variance reduction formula.
[/proofplan]
[step:Center the estimator and identify its mean]
Define the centered random variables $\widetilde Y,\widetilde C \in L^2(\Omega,\mathcal F,\mathbb P;\mathbb R)$ by
\begin{align*}
\widetilde Y := Y-\mathbb E[Y],
\end{align*}
and
\begin{align*}
\widetilde C := C-\mathbb E[C].
\end{align*}
For each $\beta \in \mathbb R$, linearity of expectation gives
\begin{align*}
\mathbb E[Y_\beta]=\mathbb E[Y]-\beta\mathbb E[C-\mathbb E[C]]=\mathbb E[Y].
\end{align*}
Therefore
\begin{align*}
Y_\beta-\mathbb E[Y_\beta]=\widetilde Y-\beta\widetilde C.
\end{align*}
Since $Y,C \in L^2(\Omega,\mathcal F,\mathbb P;\mathbb R)$ and constants belong to $L^2(\Omega,\mathcal F,\mathbb P;\mathbb R)$, the random variables $\widetilde Y$, $\widetilde C$, and $\widetilde Y-\beta\widetilde C$ are square-integrable.
[/step]
[step:Expand the variance as a quadratic in $\beta$]
For every $\beta \in \mathbb R$, the definition of variance gives
\begin{align*}
\operatorname{Var}(Y_\beta)=\mathbb E[(\widetilde Y-\beta\widetilde C)^2].
\end{align*}
Expanding the square inside the expectation and using linearity of expectation,
\begin{align*}
\operatorname{Var}(Y_\beta)=\mathbb E[\widetilde Y^2]-2\beta\mathbb E[\widetilde Y\widetilde C]+\beta^2\mathbb E[\widetilde C^2].
\end{align*}
By the definitions of variance and covariance,
\begin{align*}
\mathbb E[\widetilde Y^2]=\operatorname{Var}(Y),
\end{align*}
\begin{align*}
\mathbb E[\widetilde Y\widetilde C]=\operatorname{Cov}(Y,C),
\end{align*}
and
\begin{align*}
\mathbb E[\widetilde C^2]=\operatorname{Var}(C).
\end{align*}
Hence
\begin{align*}
\operatorname{Var}(Y_\beta)=\operatorname{Var}(Y)-2\beta\operatorname{Cov}(Y,C)+\beta^2\operatorname{Var}(C).
\end{align*}
[guided]
The goal is to make the dependence on $\beta$ explicit. Since $Y_\beta$ has the same expectation as $Y$, its centered form is especially simple:
\begin{align*}
Y_\beta-\mathbb E[Y_\beta]=\widetilde Y-\beta\widetilde C.
\end{align*}
Therefore the variance is the expectation of the square of this centered expression:
\begin{align*}
\operatorname{Var}(Y_\beta)=\mathbb E[(\widetilde Y-\beta\widetilde C)^2].
\end{align*}
Because the variables are real-valued, the usual algebraic identity $(a-b)^2=a^2-2ab+b^2$ applies pointwise, so
\begin{align*}
(\widetilde Y-\beta\widetilde C)^2=\widetilde Y^2-2\beta\widetilde Y\widetilde C+\beta^2\widetilde C^2.
\end{align*}
All three terms are integrable: $\widetilde Y^2$ and $\widetilde C^2$ are integrable by square-integrability, and $\widetilde Y\widetilde C$ is integrable by Cauchy-Schwarz in $L^2(\Omega,\mathcal F,\mathbb P)$. Applying linearity of expectation gives
\begin{align*}
\operatorname{Var}(Y_\beta)=\mathbb E[\widetilde Y^2]-2\beta\mathbb E[\widetilde Y\widetilde C]+\beta^2\mathbb E[\widetilde C^2].
\end{align*}
The definitions of variance and covariance now identify the three coefficients:
\begin{align*}
\mathbb E[\widetilde Y^2]=\operatorname{Var}(Y),
\end{align*}
\begin{align*}
\mathbb E[\widetilde Y\widetilde C]=\operatorname{Cov}(Y,C),
\end{align*}
and
\begin{align*}
\mathbb E[\widetilde C^2]=\operatorname{Var}(C).
\end{align*}
Thus the variance is the quadratic function
\begin{align*}
\operatorname{Var}(Y_\beta)=\operatorname{Var}(Y)-2\beta\operatorname{Cov}(Y,C)+\beta^2\operatorname{Var}(C).
\end{align*}
[/guided]
[/step]
[step:Complete the square to locate the unique minimizing coefficient]
Set
\begin{align*}
a := \operatorname{Var}(C),
\end{align*}
and
\begin{align*}
b := \operatorname{Cov}(Y,C).
\end{align*}
By hypothesis, $a>0$. The preceding step gives, for every $\beta \in \mathbb R$,
\begin{align*}
\operatorname{Var}(Y_\beta)=\operatorname{Var}(Y)-2b\beta+a\beta^2.
\end{align*}
Completing the square,
\begin{align*}
a\beta^2-2b\beta=a\left(\beta-\frac{b}{a}\right)^2-\frac{b^2}{a}.
\end{align*}
Therefore
\begin{align*}
\operatorname{Var}(Y_\beta)=\operatorname{Var}(Y)-\frac{b^2}{a}+a\left(\beta-\frac{b}{a}\right)^2.
\end{align*}
Since $a>0$, the last term is nonnegative for all $\beta \in \mathbb R$ and equals $0$ exactly when $\beta=b/a$. Hence the unique minimizer is
\begin{align*}
\beta_*=\frac{b}{a}=\frac{\operatorname{Cov}(Y,C)}{\operatorname{Var}(C)}.
\end{align*}
[/step]
[step:Evaluate the variance at the minimizing coefficient]
Substituting $\beta=\beta_*$ into the completed-square formula gives
\begin{align*}
\operatorname{Var}(Y_{\beta_*})=\operatorname{Var}(Y)-\frac{b^2}{a}.
\end{align*}
Using the definitions of $a$ and $b$, this is
\begin{align*}
\operatorname{Var}(Y_{\beta_*})=\operatorname{Var}(Y)-\frac{\operatorname{Cov}(Y,C)^2}{\operatorname{Var}(C)}.
\end{align*}
This proves both the optimal coefficient formula and the resulting minimum variance formula.
[/step]
