Is coupled cluster variational for two electrons?
No. The reason for this is not to be found in the excitations, but in the evaluation of the method, i.e. the working equations.
$$%Introducing some shortcuts
\require{cancel}
\newcommand{\op}[1]{\mathrm{#1}}  %\op{H}
\newcommand{\bracket}[2]{\left\langle#1\middle|#2\right\rangle}
\newcommand{\bra}[1]{\left\langle#1\right|}
\newcommand{\ket}[1]{\left|#1\right\rangle}
\newcommand{\order}[1]{^{(#1)}}         %E_n\order{1}
\newcommand{\overlap}[3]{\mathcal{#1}_{#2}\order{#3}} %\overlap{S}{m}{1}
\newcommand{\integral}[3]{\mathcal{#1}_{#2,#3}} %\integral{V}{i}{j}
$$

The variational principle

In short: a trial wave function always has a higher energy than the exact wave function.

Assume we know the exact solutions to the Schrödinger equation, which form a complete basis and there is one wave function $\Psi_0$ which leads to the lowest energy $E_0$.
\begin{align}
\op{H}\Psi_i &amp;= E_i\Psi_i &amp; 
\{i \in \mathbb{N}, i\geq 0\}\tag{1}\label{schrodinger}
\end{align}

We choose the solutions to be orthonormal (boundary conditions).
$$\bracket{\Psi_i}{\Psi_j} = \delta_{ij} = 
\begin{cases}1, &amp; i=j \\ 0, &amp; i\neq j \end{cases}
\tag2\label{orthonorm}$$

The trial wave function $\Phi$ can be expressed as a linear combination of the complete set of wave functions.
$$\Phi = \sum_i^\infty a_i\Psi_i\tag3\label{expansion}$$

The expectation energy of the trial wave function is given through
$$W = \frac{\bra{\Phi}\op{H}\ket{\Phi}}{\bracket{\Phi}{\Phi}}.
\tag4\label{expen}$$

With $\eqref{expansion}$, $\eqref{orthonorm}$, and $\eqref{schrodinger}$ we can expand and simplify.
\begin{align}
W &amp;=
  \frac{
    \sum_i^\infty\sum_j^\infty a_i a_j \bra{\Psi_i}\op{H}\ket{\Psi_j}
  }{
    \sum_i^\infty\sum_j^\infty a_i a_j \bracket{\Psi_i}{\Psi_j}
  }\\
  &amp;=
  \frac{
    \sum_i^\infty\sum_j^\infty a_i a_j E_j\bracket{\Psi_i}{\Psi_j}
  }{
    \sum_i^\infty\sum_j^\infty a_i a_j \bracket{\Psi_i}{\Psi_j}
  }\\
  &amp;=
  \frac{
    \sum_i^\infty a_i^2 E_i\bracket{\Psi_i}{\Psi_j} +
    \sum_{j\neq i}^\infty a_i a_j E_j\bracket{\Psi_i}{\Psi_j}
  }{
    \sum_i^\infty a_i^2 \bracket{\Psi_i}{\Psi_j} +
    \sum_{j\neq i}^\infty a_i a_j \bracket{\Psi_i}{\Psi_j}
  }\\
  &amp;=
  \frac{
    \sum_i^\infty a_i^2 E_i
  }{ 
    \sum_i^\infty a_i^2
  } &amp; | -E_0\\
W - E_0 &amp;=
  \frac{
    \sum_i^\infty a_i^2 (E_i - E_0)
  }{ 
    \sum_i^\infty a_i^2
  } \\
W - E_0 &amp;=
  \frac{
    \sum_{i&gt;0}^\infty a_i^2 (E_i - E_0)
  }{ 
    \sum_{i&gt;0}^\infty a_i^2
  } + \frac{a_0^2 (E_0 - E_0)}{a_0^2}\tag{4'}\label{notnegative}\\
W - E_0 &amp;\geq 0\tag5\label{varprin}
\end{align}

In $\eqref{notnegative}$ the last term is zero. Since $E_i &gt; E_0$ (by definition) and $a_i^2$ is always greater or equal zero, the expectation value of the trial wave function must always have a higher energy than the exact solution of the Schrödinger equation, i.e. $\eqref{varprin}$.

Configuration Interaction

The CI wave function is set up as a linear combination of determinants from a reference calculation. This reference calculation can in principle be any other method the only requirement is that it also forms a complete (orthonormal) basis.
When Hartree-Fock is chosen, we know that we obtain a complete set and that this method itself is variational. The CI trial wave function can then be expressed as a linear combination of solutions of the HF wave function.
$$\Psi^\mathrm{CI} = \sum_i^\infty a_i\Phi_i = 
a_0\Phi^\mathrm{HF} + \sum_{i&gt;0}^\infty a_i\Phi_i$$
An alternative expression of the above is using an "excitation operator" $\op{T}$ to generate the higher energy determinants from the reference.
$$\Psi^\mathrm{CI} = (1 + \op{T})\Phi^\mathrm{HF}$$

The energy of the CI trial wave function can be evaluated with the Lagrange multiplier method, which requires, that the wave function is normalised (adding zero).
$$L = \bra{\Psi^\mathrm{CI}}\op{H}\ket{\Psi^\mathrm{CI}} 
     - \lambda\left[\bracket{\Psi^\mathrm{CI}}{\Psi^\mathrm{CI}}-1\right]$$

Since the boundary conditions did not change, the variational principle still holds.

Coupled Cluster

Because of the inherent complexity of the CC ansatz, one key element to evaluate to CC energy is missing: the requirement of the normalised CC trial wave function.
We generate the trial wave function with an exponential approach (rather than a linear) expanded into a Taylor series.
\begin{align}
\Psi^\mathrm{CC} 
&amp;= \mathrm{e}^{\op{T}}\Phi^\mathrm{HF}\\
&amp;= \left[\op{1} + \op{T}_1 + \left(\op{T}_2 + \frac12\op{T}_1^2\right)
   + \left(\op{T}_3 + \op{T}_2\op{T}_1 + \frac16\op{T}_1^3\right)
   + \cdots \right] \Phi^\mathrm{HF}\\
&amp;= \Phi^\mathrm{HF} + \left[\op{T}_1 + \left(\op{T}_2 + \frac12\op{T}_1^2\right)
   + \left(\op{T}_3 + \op{T}_2\op{T}_1 + \frac16\op{T}_1^3\right)
   + \cdots \right] \Phi^\mathrm{HF}\\
\end{align}

Or to shorten the journey
$$\Psi^\mathrm{CC} = \Phi^\mathrm{HF} + \sum_i^\infty t_i \Phi_i.$$

One can easily see what kind of mess it would become requiring the wave function to be normalised $\bracket{\Psi^\mathrm{CC}}{\Psi^\mathrm{CC}}=1$, so the Lagrange ansatz won't work. Instead we choose the wave function to be orthonormal to the reference, which in the Hartree-Fock (complete orthonormal basis) takes away a lot of work.
\begin{align}
\bracket{\Phi^\mathrm{HF}}{\Psi^\mathrm{CC}}&amp;= 1&amp;
\implies\quad \bracket{\Phi^\mathrm{HF}}{\Phi_i} &amp;= 0
\end{align}

We can write the CC-Schrödinger equation as
$$\op{H}\mathrm{e}^{\op{T}}\Phi^\mathrm{HF} = 
  E^\mathrm{CC}\mathrm{e}^{\op{T}}\Phi^\mathrm{HF}$$
and evaluate the energy as 
\begin{align}
\bra{\Phi^\mathrm{HF}}\op{H}\mathrm{e}^{\op{T}}\ket{\Phi^\mathrm{HF}}
  &amp;= E^\mathrm{CC}
     \bracket{\Phi^\mathrm{HF}}{\mathrm{e}^{\op{T}}\Phi^\mathrm{HF}}\\
E^\mathrm{CC} 
  &amp;= \bra{\Phi^\mathrm{HF}}\op{H}\mathrm{e}^{\op{T}}\ket{\Phi^\mathrm{HF}}\\
E^\mathrm{CC} 
  &amp;= E^\mathrm{HF} + 
  \bra{\Phi^\mathrm{HF}}\op{H}
    \ket{\left[\op{T}_1 + \left(\op{T}_2 + 
     \frac12\op{T}_1^2\right)
     + \left(\op{T}_3 + \op{T}_2\op{T}_1 + \frac16\op{T}_1^3\right)
     + \cdots \right]\Phi^\mathrm{HF}}.\\
\end{align}

From this we can see that the evaluation of the expectation value of the energy is significantly different from what is required for the variational principle $\eqref{expen}$.

Without going further into details, the dropping of the requirement of a normalised trial wave function leads to not obeying the variational principle.

Additional considerations

It is often said that Full-CC is equivalent to Full-CI, which is a convenient half truth. It is extremely important that this gets more true with a better description of the reference wave function, i.e. with larger basis sets. Full-CI is exact at the complete basis set limit (for the electronic Schrödinger equation). This is not necessarily true for Full-CC. 

Since we are in computational chemistry hardly ever concerned about absolute energies, the non-variational nature of the method is by far outperformed by the "infinite" approach generating the "excited" determinants. Relative energies don't suffer significantly (compared with the size consistency and extensivity) under the non-variational approach.