How do I correctly interpret $\rho = \psi_1^* \psi_2$? 
Summary: This turned out to be a rather trivial one indeed. As Marek mentioned in the comment, the continuity equation is trivial. And it indeed turns out be so. Godfrey Miller elaborates on this, showing that the continuity equation merely shows that the transition amplitude is constant with time. All this is only confusing as long as you stick to the belief that  the continuity equation is some exclusive indication of probability density. 

$\rho = \psi_1^* \psi_2$ is the probability density for transition from state 2 to 1.  I am having trouble interpreting the probability though $$dP_{1\leftarrow 2}=\psi_1^* \psi_2 d^3 r $$ does this mean the probability in a region of volume $d^3 r$ to go from state 2 to state 1? That does not make sense at all!
EDIT: I more or less agree with Marek's answer. However, that makes it seem this is not the probability density. How do we reconcile that with the fact that it satisfies the continuity equation. If we have two states:
$$ i \hbar \partial_t \psi_{1,2}= -\frac{\hbar^2}{2m} \nabla^2 \psi_{1,2} + V\psi_{1,2}$$
After making standard manipulations (taking conjugate, multiplying by conjugate, subtracting) we would get
$$\partial_t ( \psi_1^*\psi_2) +\frac{\hbar}{2mi} \nabla\cdot \left( \psi_1^* \nabla \psi_2 - (\nabla \psi_1^* )\psi_2\right)=0$$
This seems consistent with the continuity equation: $$\partial_t \rho +\nabla.\mathbb{j} =0$$ However, we seem to be changing its relation to the probability by making the probability $P=\int \int\rho^*(r,r') \rho(r,r') d^3r d^3r' $ instead of the regular $P=\int \rho(r) d^3r$ For a kind of intuitive "dimensional consistency" can we instead say $P=\int \int \; \sqrt{\rho^* \rho} \; d^3r d^3r'$
$$****$$
Full steps: 
$\psi_1$ and $\psi_2$ are two states of the same system, i.e both satisfy the SE of the same potential:
Taking the SE for $\psi_2$  and multiplying by $\psi_1^*$ :
$$ \psi_1^*\times\left( i \hbar \partial_t \psi_2 = -\frac{\hbar^2}{2m} \nabla^2 \psi_2 + V\psi_2\right)$$ Taking the SE for $\psi_1$ complex conjugating the entire equation and multiplying by $\psi_2$ :
$$\left( -i \hbar \partial_t \psi_1^*= -\frac{\hbar^2}{2m} \nabla^2 \psi_1^* + V\psi_1^*\right)\times \psi_2$$
($V$ is real, no decay processes etc) On subtracting the second equation from the first 
$$\partial_t (\psi_1^*\psi_2 )= -\frac{\hbar}{2mi} \left( \psi_1^* \nabla^2 \psi_2 - (\nabla^2 \psi_1^* )\psi_2\right)$$
 which is the same as $$ \partial_t (\psi_1^*\psi_2 )= -\frac{\hbar}{2mi} \nabla\cdot \left( \psi_1^* \nabla \psi_2 - (\nabla \psi_1^* )\psi_2\right)$$
 A: That's probably (no pun intended) because it's not a probability density at all. Wave functions by themselves don't have any physical interpretation. Only their squares do.
The probability of transition is computed like this
$$P_{1 \to 2} = {\left \Vert \left < \psi_1 | \psi_2 \right > \right \Vert^2 \over \left \Vert \psi_1 \right \Vert^2 \left \Vert \psi_2 \right \Vert^2}$$
What you've written is just a density of the probability amplitude $\left<\psi_1|x\right>\left<x|\psi_2\right>$. It can be complex so it certainly cannot be probability in any sense.
Edit regarding your additional comments: assuming that $\left \Vert \psi_1 \right \Vert^2 = \left \Vert \psi_2 \right \Vert^2 = 1$ you could write $P_{1 \to 2}$ like this 
$$P_{1 \to 2} = {\rm Tr} \left( \left|\psi_1\left> \right<\psi_1\right|  \left|\psi_2\left> \right<\psi_2\right| \right) = \int \psi_1(y) \psi_1^*(x)  \psi_2(x) \psi_2^*(y) {\rm d} x {\rm d} y$$
interpreting it as a trace of the product of projection operators. But this integral representation is basically never done, it's much easier to only compute the amplitude and square it when you're finished.
As a moral, you shouldn't insist on restricting your attention to some portion of position space. Quantum theory "doesn't like it" (except for few special constructions).
A: Since
\begin{equation}
\langle \psi_1 \mid \psi_2 \rangle = \int d^3 x \ \psi_1^*(t,x) \psi_2(t,x),
\end{equation}
the quantity $\psi_1^* \psi_2$ is the spatial density of the transition amplitude from $\psi_2$ to $\psi_1$.  This amplitude density is in general complex, so it is not possible in general to interpret it as a probability density, but according to your equation it is conserved all the same.  The only exception is when $\mid \psi_1 \rangle = \mid \psi_2 \rangle = \mid \psi \rangle$, in which case the equation above becomes
\begin{equation}
1 = \int d^3 x \ \psi^*(t,x) \psi(t,x),
\end{equation}
which is why the non-negative quantity $\psi^* \psi$ can be interpreted as a probability density. 
The conservation of amplitude density implies the important fact that transition amplitudes between well-behaved solutions of the Schrodinger equation are time independent.  It follows from your equation that
\begin{equation}
\frac{\partial}{\partial t} \langle \psi_1 \mid \psi_2 \rangle = \int dx \ \frac{\partial}{\partial t} (\psi_1^* \psi_2) = \frac{i \hbar}{2 m} \int dx \ \nabla \cdot (\psi_1^* \nabla \psi_2 - \psi_2 \nabla \psi_1^*).
\end{equation}
If the wavefunctions $\psi_1$, $\psi_2$ are well-behaved at infinity, this final integral will vanish, so $\langle \psi_1 \mid \psi_2 \rangle = {\rm constant}$.  In hindsight, you could have inferred the necessity of the conservation of amplitude density from two facts: 1) transition amplitudes can be written as an integral of a local quantity, and 2) transition amplitudes are time-independent.
Does this help answer your question?
A: I'm not sure what a "probability density for transition" is - is it a probability density or a transition amplitude?  I think you mean transition amplitude.
Are you sure that $\rho = \psi_1^* \psi_2$ is a transition amplitude?  It looks more like simply an overlap.  How can it be a transition amplitude without any time evolution operator or Hamiltonian?  I believe that a transition amplitude should have an interaction which connects the two states appearing between them.
If I assume that $\psi_1$ and $\psi_2$ separately satisfy Schrodinger's Equation, then the quantity $\psi_1^* \psi_2$ would be a conserved charge with the current calculated above. This "conservation law" is just a statement that the overlap between two states satisfying Schrodinger's Equation will be preserved as the states evolve in time.
If I assume that $\psi_1$ and $\psi_2$ separately satisfy Schrodinger's Equation and the matrix element between the states is 1, then
\begin{equation}
\int d^3 x \ \psi_1^*(x) \psi_2(x)
\end{equation}  
is a transition amplitude, I think.  The probability of transition would then go like the modulus of the amplitude squared.  This is different from a probability density.
Are you conflating two ideas?  Have I missed something?
A: Even the scalar product $\langle\psi_1|\psi_2\rangle$ is not a transition probability but an amplitude to find the system in one state when it is know that it is in another. And the wave function product is just a part of this amplitude. It's like a single term of a scalar product of two vectors $(\vec{a}\vec{b})$: $a_i b_i$, $i$ = fixed.
EDIT: "Transition" implies evolution in time. It is not the case here. It is amplitude of "presence". Although the wave functions may evolve in time (each is a non trivial superpositions of eigenstates with their "energetic" exponentials) the time evolution of each state is determined with the same Hamiltonian, so the product $\psi_1 (t)^*\psi_2 (t) = \psi_1 (0)^*\psi_2 (0)$ due to $U(t)^+U(t)=1$.
