Maybe this little "proof" does clear some things for you. I especially think that you do not quite understand the correlation coefficient entirely, but thats just my guess and I don't want to offend you.
First let me remind you that the correlation coefficient is a measure of the linear dependency of two independent random variables X and Y.
So lets look at the formula for the covariance and the correlation coefficient:
Cov(X,Y) = E[(X-E[X])(Y-E[Y])
$r_{XY} = \frac{Cov(X,Y)}{\sqrt{Var(X)}\sqrt{Var(Y)}}$
Since X,Y ~ N(0,1) we get: $Cov(X,Y) = r_{XY}$
So now lets get to linear regression:
$Y = \beta_{0} + \beta_{1}X + \epsilon$
Notice that the $\beta_{0}$ intercept is cancelled out, because of X,Y ~ N(0,1), we get:
$Y = \beta_{1}X + \epsilon$
Now note that the slope $\beta_{1}$ is defined as: $r_{XY} * \frac{\sqrt{Var(Y)}}{\sqrt{Var(X)}}$, since X,Y ~ N(0,1):
$\beta_{1} = r_{XY}$ = Cov(X,Y)
So now lets move on to your example:
What you want to fix is:
$A = r_{AB}B + \epsilon_{AB}$
$P = r_{BP}B + \epsilon_{BP}$
What you want to "estimate" is
$r_{AP}$, which is equal to Cov(A,P) since A,P ~ N(0,1) as shown before.
The Covariance formula can be rewritten as:
Cov(X,Y) = E[XY] - E[X]*E[Y]
So we get $E[AP] - \underbrace{E[A]}_{0}\underbrace{E[B]}_{0}$, again since A,B ~ N(0,1)
Cov(A,P) = E[AP]
Lets plug our two equations for A and P in there:
$E[AP] = E[(r_{AB}B + \epsilon_{AB})(r_{BP}B + \epsilon_{BP}])]$
$E[AP] = E[r_{AB}B^{2}r_{BP}] + \underbrace{E[r_{AB}B\epsilon_{BP}]}_{0} + \underbrace{E[\epsilon_{AB}Br_{BP}]}_{0} + \underbrace{E[\epsilon_{AB}\epsilon_{BP}]}_{0}$, since we assume uncorrelated residuals $\epsilon$ (crucial regression assumption)
$E[AP] = r_{AB}r_{BP}E[B^{2}]$, because $E[r_{AB}] = r_{AB}$ and $E[r_{BP}] = r_{BP}$
$E[AP] = r_{AB}r_{BP}\underbrace{E[B^{2}]}_{1}$, because B ~ N(0,1)
so we finally get $E[AP] = r_{AB}r_{BP}$, which can be rewritten as:
$r_{AP} = r_{AB}r_{BP}$
So now you also should get an idea of when this not holds true. For example if the residuals $\epsilon$ are not uncorrelated (In reality they are almost never entirely uncorrelated).
Maybe @Glen_b♦ can check this answer for mistakes or advance further on situations, when this does not holds true, because he has far more knowledge than I have.
Suppose $X$ is multivariate normal.
Then by definition we can write $X = \mu + CZ$ for a (deterministic) vector $\mu$, matrix $C$ and random vector $Z$ of independent standard normal (i.e. $N(0,1)$) variables.
(So $\Sigma = CC^T$ is the covariance matrix).
Thus if $C$ is invertible then
$$
C^{-1}(X-\mu) = Z
$$
is a vector of independent (and hence uncorrelated) standard normals.
Note, if you are given pos. def. $\Sigma$ upfront, then to calculate $C = (\Sigma^{1/2})^{-1}$ you would probably want to first calculate the square root (e.g. using Cholesky) and then invert that.
Best Answer
Binomial variables are usually created by summing independent Bernoulli variables. Let's see whether we can start with a pair of correlated Bernoulli variables $(X,Y)$ and do the same thing.
Suppose $X$ is a Bernoulli$(p)$ variable (that is, $\Pr(X=1)=p$ and $\Pr(X=0)=1-p$) and $Y$ is a Bernoulli$(q)$ variable. To pin down their joint distribution we need to specify all four combinations of outcomes. Writing $$\Pr((X,Y)=(0,0))=a,$$ we can readily figure out the rest from the axioms of probability: $$\Pr((X,Y)=(1,0))=1-q-a, \\\Pr((X,Y)=(0,1))=1-p-a, \\\Pr((X,Y)=(1,1))=a+p+q-1.$$
Plugging this into the formula for the correlation coefficient $\rho$ and solving gives $$a = (1-p)(1-q) + \rho\sqrt{{pq}{(1-p)(1-q)}}.\tag{1}$$
Provided all four probabilities are non-negative, this will give a valid joint distribution--and this solution parameterizes all bivariate Bernoulli distributions. (When $p=q$, there is a solution for all mathematically meaningful correlations between $-1$ and $1$.) When we sum $n$ of these variables, the correlation remains the same--but now the marginal distributions are Binomial$(n,p)$ and Binomial$(n,q)$, as desired.
Example
Let $n=10$, $p=1/3$, $q=3/4$, and we would like the correlation to be $\rho=-4/5$. The solution to $(1)$ is $a=0.00336735$ (and the other probabilities are around $0.247$, $0.663$, and $0.087$). Here is a plot of $1000$ realizations from the joint distribution:
The red lines indicate the means of the sample and the dotted line is the regression line. They are all close to their intended values. The points have been randomly jittered in this image to resolve the overlaps: after all, Binomial distributions only produce integral values, so there will be a great amount of overplotting.
One way to generate these variables is to sample $n$ times from $\{1,2,3,4\}$ with the chosen probabilities and then convert each $1$ into $(0,0)$, each $2$ into $(1,0)$, each $3$ into $(0,1)$, and each $4$ into $(1,1)$. Sum the results (as vectors) to obtain one realization of $(X,Y)$.
Code
Here is an
Rimplementation.