I think the simplest way to understand why we cannot have a joint probability distribution for two incompatible observables $A,B$ (meaning non-commuting) is the following. With a slight abuse of of notation, the joint probability is defined as:

$$
P_{A,B}(a,b) := \mathrm{Prob}(A=a,B=b)
$$

which means, it's the probability of $A$ having value $a$ and $B$ having value $b$. In quantum mechanics this means that there is an eigenstate of $A$ with eigenvalue $a$ and an eigenstate of $B$ with eigenvalue $b$. But if $[A,B]\neq 0$ this is notoriously not possible. 

Note that, conversely, if $A$ and $B$ commute it is always possible to find a common eigenbasis and so the above prescription works fine. 

This argument, however does not answer the other part of the question. Which is:
> Why can I not define a *bona fide* characteristic function (i.e. which is the Fourier transform of a probability density) in case of non-compatible observables? And perhaps, is there any way to amend this? 

As pointed correctly by the OP this has to do with Bochner's theorem which tells precisely what requirements a characteristic function has to satisfy. I will state Bochner's theorem in the form needed for our purposes.

> **Theorem (Bochner)** $\chi (t)$ is the Fourier transform of a probability density $P(\omega)$ if and only if for any $n$-tuple $t_1,t_2,\ldots t_n$ the matrix with entries $\chi_{i,j}:= \chi(t_i-t_j)$ is non-negative definite (and hermitian).

A simple way to understand Bochner's theorem is the following. A matrix $\chi$ is non-negative definite if and only if it can be written as $\chi = A A^\dagger$.

Let $P(\omega)$ be the Fourier transform of $\chi$, that is 

$$
\chi_{i,j} := \int d\omega e^{i(t_i-t_j) \omega} P(\omega)
$$

which we write as

\begin{align}
\chi_{i,j} &:= \int d\omega e^{it_i \omega} P(\omega) e^{-it_j \omega} \\
& = (A A^\dagger)_{i,j}
\end{align}

with

$$
A_{i,\omega} := e^{it_i \omega} \sqrt{P(\omega)}
$$

which we can do since $P(\omega)$ is non-negative. So $P(\omega)$ non-negative means that $\chi_{i,j}$ is a non-negative definite matrix. This characterizes characteristic functions in the classical case. 

Let us now turn to the quantum mechanics and consider the multivariate case, i.e., we have several observables which I call $X_1, X_2, \ldots X_n$ with spectra in $\omega_1, \ldots, \omega_n$. The conjugate variables being $t_1,\ldots,t_n$, and the notation 

$$tX:=\sum_{k=1}^n t_k X_k$$

In this case the characteristic function is

$$
\chi(t):=\mathsf{E}[ e^{itX} ]= \operatorname{Tr} ( e^{itX} \rho )
$$

for some quantum state $\rho$ (a normalized non-negative matrix). We want to check under which conditions $\chi_{i,j}=\chi(t_1-t_j)$ is non-negative definite as a matrix. 

If the $X_k$ were mutually commuting operators we would have

$$
e^{i(t_i-t_j) X} = e^{it_i X} e^{-it_j X} \ \ \ \ \ \ (1)
$$

and then we could write

\begin{align}
\chi_{i,j} &=\operatorname{Tr} \left ( e^{i(t_i-t_j) X} \rho \right)  \\
&=\operatorname{Tr} \left ( e^{it_i X} e^{-it_j X}  \rho \right )  \\
&=\operatorname{Tr}  \left  (e^{-it_j X}  \sqrt{\rho} \sqrt{\rho} e^{it_i X} \right )  
\end{align}

Now define the matrix $A_{j,lq} := \left ( e^{-it_j X}  \sqrt{\rho} \right )_{l,q} =: A_{i,\xi}$ where $\xi=(l,q)$. We have

$$
\chi_{i,j} = \sum_{lq} A_{j,lq} \overline{A_{j,lq}}
$$

which is of the form $BB^\dagger$ and proves that $\chi_{i,j}$ is non-negative definite. 
Obviously this whole construction breaks down if $X_k$ are not mutually commuting. Since the condition in Bochner's theorem is an if and only if, in this case $\chi(t)$ is not the Fourier transform of a probability distribution