Proof of $P(\mathbf a,\mathbf b)=-\mathbf a\cdot\mathbf b$ in Bell’s original paper Bell’s original 1964 argument works by comparing a QM prediction with the prediction of a local hidden variable theory. In the notation of Griffith’s version of it in ”Introduction to Quantum Mechanics” where $P(\mathbf{a},\mathbf{b})$ is the average value of the product of the spins for a given set of detector orientations given by the unit vectors $\mathbf{a},\mathbf{b}$, the QM prediction is:
$P(\mathbf{a},\mathbf{b}) = -\mathbf{a} \cdot \mathbf{b}$.
Is there a simple proof of this prediction, suitable for a student at the level of Griffith’s book? 
 A: One route is the one followed in the notes already linked in the comments. Here is another derivation (though most likely not the one used in Bell's original paper).
We want to compute the expectation value:
$$
  \langle \Psi|
    [(\boldsymbol\sigma\cdot\mathbf a) \otimes
    (\boldsymbol\sigma\cdot\mathbf b)]
  |\Psi\rangle=
  a_i b_j \Psi^\dagger (\sigma_i\otimes\sigma_j)\Psi =
  a_i b_j \bar\Psi_{\alpha\beta} \Psi_{\gamma\delta} (\sigma_i)_{\alpha\gamma}(\sigma_j)_{\beta\delta}.
$$
Now, if $|\Psi\rangle$ is maximally entangled, and in particular is the singlet $|01\rangle-|10\rangle$, then $\Psi_{\alpha\beta}=(-1)^{\alpha}\delta_{\alpha+\beta=0}/\sqrt 2=(iY)_{\alpha\beta}/\sqrt2$, where $Y\equiv\sigma_2$ is the Pauli matrix, and thus
$$
  a_i b_j \bar\Psi_{\alpha\beta} \Psi_{\gamma\delta} (\sigma_i)_{\alpha\gamma}(\sigma_j)_{\beta\delta} =
  \frac12 a_i b_j \operatorname{Tr}(Y\sigma_i Y\sigma_j^T) 
  = -\frac12 a_i b_j  \operatorname{Tr}(\sigma_i\sigma_j)
  = -\mathbf a\cdot\mathbf b,
$$
where we used the identities:
$$Y\sigma_i Y=(-1)^{1+\delta_{i,2}}\sigma_i,
\qquad \sigma_i^T = (-1)^{\delta_{i,2}} \sigma_i.$$

More generally, suppose $|\Psi\rangle$ is maximally entangled. That means that $\operatorname{Tr}_B(|\Psi\rangle\!\langle\Psi|)=I/2$, that is, in terms of the matrix elements $\Psi_{12}$, $\Psi_{12}\bar\Psi_{32}=\delta_{13}/2$, or equivalently $\Psi\Psi^\dagger =I/2$.
The general form of $\Psi$ is therefore
$$\sqrt2 \,|\Psi\rangle= |0\rangle\!\langle u| + |1\rangle\!\langle v|,$$
for some pair of orthonormal states $|u\rangle,|v\rangle$, and thus in matrix notation, $\sqrt2\,\Psi=U$, with $U$ the unitary whose rows are (the complex conjugates of) $|u\rangle$ and $|v\rangle$.
Now, any $2\times2$ unitary can be written in terms of Pauli matrices as
$$U=e^{i\theta} (c_0 I + ic_1 X +ic_2 Y + ic_3 Z),$$
for some choice of real parameters $c_0,c_1,c_2,c_3$ that are normalised to one: $(c_0,c_1,c_2,c_3)\in S^3$, and $\theta\in\mathbb R$. We can fix $\theta=0$ without loss of generality here.
Our expectation value then reads
$$
  \langle \Psi|
    [(\boldsymbol\sigma\cdot\mathbf a) \otimes
    (\boldsymbol\sigma\cdot\mathbf b)]
  |\Psi\rangle=
  a_i b_j \bar\Psi_{\alpha\beta} \Psi_{\gamma\delta} (\sigma_i)_{\alpha\gamma}(\sigma_j)_{\beta\delta} \\
  = a_i b_j \operatorname{Tr}(\Psi^\dagger\sigma_i \Psi \sigma_j^T)
  = \frac12 a_i b_j \operatorname{Tr}(U^\dagger\sigma_i U \sigma_j^T).
$$
There is then always an orthogonal transformation $B_{ij}$ such that $U^\dagger \sigma_i U=B_{ji}\sigma_i$, using which we conclude that
$$\frac12 a_i b_j B_{ji} \operatorname{Tr}(\sigma_i \sigma_j^T)
= \frac12 \mathbf b \tilde B \mathbf a,$$
where $\tilde B$ differs from $B$ only in the sign of its second row.
A: Let me make an attempt at a high-level argument:
We want to compute 
$$
\langle\Psi^-\vert(\vec a \cdot\vec\sigma)\otimes(\vec b\cdot\vec \sigma)\vert\Psi^-\rangle\ .\tag{*}
$$
with $\vert\Psi^-\rangle$ the singlet state.
First, since $\vert\Psi^-\rangle$ is invariant under rotations, (*) can only depend on the relative angle between $\vec a$ and $\vec b$, that is, $\vec a \cdot \vec b$.  In addition, (*) is linear in $\vec a$ and $\vec b$.  Thus, it must be proportional to $\vec a\cdot \vec b$. The proportionality constant can now be determined by picking a specific case (such as parallel $\vec a$ and $\vec b$).
