Fortunately I have received some explanations on these questions from Slava Rychkov. So here is what I understood from his arguments on unitary CFT, plus speculations of my own on non-unitary CFT.

**A general proof**

First of all, some notations. Consider a QFT with a space of states $S$, together with a sesquilinear form (i.e. a scalar product) which we do not yet assume to be positive definite. Let us denote states as $|i\rangle \in S$ and the scalar product as $\langle i| j\rangle$. The Hermitian conjugate of an operator $O$ is by definition the operator $O^\dagger$ such that

\[

\langle i| O^\dagger |k\rangle = \langle k | O |i\rangle^*\ ,

\]

where the star denotes complex conjugation. In particular the identity is self-conjugate i.e. Hermitian. Finally we denote the state-operator correspondence as $|i\rangle \mapsto O_i$, and the Hermitian conjugate of such an operator as $O_i^\dagger = O_{i^\dagger}$.

If in addition our QFT is a CFT, then three-point functions are determined up to scalar factors by conformal symmetry. These scalar factors are called three-point structure constants, and denoted

\[

C_{ijk} = \langle O_i O_j O_k \rangle\ .

\]

We want to know how these structure constants behave under Hermitian conjugation of the operators. To do this, we use the state-operator correspondence and rewrite

\[

\langle O_i^\dagger O_j O_k \rangle = \langle i|O_j|k\rangle\ .

\]

This identity is particularly obvious if we assume that there exists a vacuum state such that $|k\rangle = O_k|0\rangle$ and $\langle \cdots \rangle = \langle 0|\cdots |0\rangle$, but this assumption is not really necessary.

From the three equations displayed so far, and from the symmetry of $C_{ijk}$ under permutations, we deduce

\[

C^*_{ijk} = C_{i^\dagger j^\dagger k^\dagger}\ ,

\]

and in particular a structure constant involving three Hermitian operators must be real.

So, if we want real structure constants, it is enough to consider Hermitian operators. But is there a Hermitian basis of operators? The answer is yes: given an operator $O$, either $O^\dagger = \lambda O$ in which case a rescaling makes it Hermitian, or $O$ and $O^\dagger$ are linearly independent in which case they can be replaced with the Hermitian combinations $O+O^\dagger$ and $i(O-O^\dagger)$.

**The result seems too strong! Where's the mistake?**

At this point the reader may protest that the argument does not rely on unitarity at all, and we seem to have proved that three-point structure constants are real in any CFT. This is actually true, but not as strong as it appears. The point is that Hermitian operators need not have positive two-point functions, except in unitary CFTs. Given a Hermitian operator $O_i$, we indeed have

\[

\langle O_i O_i \rangle = \langle O_i^\dagger O_i \rangle = \langle i | i \rangle \ ,

\]

which is always real, and positive in a unitary CFT. In a non-unitary CFT, insisting on the fundamental assumption that two-point functions are positive (or even normalized to one) can make structure constants imaginary. So we can tentatively formulate the result as

In any CFT, there are Hermitian bases of operators, and the corresponding three-point structure constants are real. In a unitary CFT, moreover, the two-point functions are positive.

**The case of Minimal Models**

Two-dimensional CFT provides a wealth of fully solved CFTs, unitary or not. So this is a good place to look for (counter)examples. Let us consider the A-series Virasoro Minimal Models (reviewed here). The simplest non-unitary model is the model of the Lee-Yang singularity $M_{5,2}$. In this model there is only one primary field $\phi$ which is not the identity, and the structure constant $C_{\phi\phi\phi}$ is purely imaginary. This suggests that this field is anti-Hermitian. Considering the Hermitian field $i\phi$ instead, we would have real structure constants, but a negative two-point function $\langle (i\phi)(i\phi) \rangle$.

More generally, whenever three-point functions are not real in Minimal Models, this is because some primary fields are anti-Hermitian. These fields can be renormalized and made Hermitian, but their two-point functions then become negative. So the example of Minimal Models supports our general claims.

Interestingly, non-unitary Minimal Models always have primary fields with negative conformal dimensions. The presence of negative dimensions is a sufficient condition for non-unitarity, and I wonder under which assumptions it is also necessary. (It cannot be necessary in general, cf the $H_3^+$ model.)

**Numerical bootstrap for non-unitary theories?**

Adapting the numerical bootstrap to non-unitary theories may not be very difficult. We can work with Hermitian operators, in order to have real structure constants. In particular, we can assume all primary operators to be Hermitian. The question is whether the two-point functions are positive or negative. There, a binary choice has to be made for each primary operator. (The choice can become more complicated if there are several primary operators with the same conformal dimension.) The case of Minimal Models suggests that the space of non-unitary CFTs may not be too much larger than the space of unitary CFTs.