I am looking for the solution to:
\begin{equation} E[VP|V^{*}<a]=?? \end{equation}
\begin{equation} \begin{pmatrix} P \\ V^{*} \end{pmatrix}\sim N\left(\begin{pmatrix} \mu_{p} \\ \mu_{v^{*}} \end{pmatrix},\begin{pmatrix} \sigma_{p} & \sigma_{v^{*},p} \\ \sigma_{v^{*},p} & \sigma_{v^{*}} \end{pmatrix}\right) \end{equation}
Where V* and P are bivariately normally distributed (do not assume independence). And
\begin{equation}V=D_{v>a}V^{*}+(1-D_{v>a})a \end{equation}
With D=1 if the condition is met, and 0 otherwise. We can interpret this as: whenever v is lower than a, it becomes a instead. This is a censored, not truncated, process.
My attempt at finding a solution is as follows, start with the general formula for the bivariate unconditional expectation for the latent variable:
\begin{equation} E[V^{*}P]=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} V^{*}P \phi{(v^{*},p)} dp dv^{*} \end{equation}
Introduce substitution, with the indicator dummy.
\begin{equation} E[VP]=E[(D_{v>a}V^{*}+(1-D_{v>a})a)P]=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} (D_{v>a}V^{*}+(1-D_{v>a})a)P \phi{(v^{*},p)} dp dv^{*} \end{equation}
Since (thx Daniel)!:
\begin{equation}E[VP \mid V^{\ast} < a] = \frac{E[VP \cdot\mathbb{1}_{V^{\ast} < a}]}{P[V^{\ast} < a]} \end{equation}
If we condition, we can rewrite as:
\begin{equation} E[VP|V^{*}<a]=\frac{1}{P[V^{\ast} < a]} \int_{-\infty}^{a} \int_{-\infty}^{\infty} (D_{v^{*}>a}V^{*}+(1-D_{v^{*}>a})a)P \phi{(v^{*},p)} dp dv^{*} \end{equation}
Since under the integration bounds, v* is always smaller than a, we can simplify:
\begin{equation} E[VP|V^{*}<a]=\frac{1}{P[V^{\ast} < a]}\int_{-\infty}^{a} \int_{-\infty}^{\infty} aP \phi{(v^{*},p)} dp dv^{*} \end{equation}
Using, the fact that a is a constant, and the fact that the joint is the product of marginal and the univariate.
\begin{equation} E[VP|V^{*}<a]=\frac{a}{P[V^{\ast} < a]} \int_{-\infty}^{a} \phi(v^{*}) \int_{-\infty}^{\infty} P \phi{(p|v^{*})} dp dv^{*} \end{equation}
Here the right term is the expectation of the marginal, which we can substitute.
\begin{equation} E[VP|V^{*}<a]=\frac{a}{P[V^{\ast} < a]} \int_{-\infty}^{a} \phi(v^{*}) [\mu_{p} + \frac{\sigma_{v^{*},p}}{\sigma_{v^{*}}^2}*(v^{*} -\mu_{v^{*}})] dv^{*} \end{equation}
Re-arranging, using linearity of integrals:
\begin{equation} E[VP|V^{*}<a]=\frac{a}{P[V^{\ast} < a]} \int_{-\infty}^{a} \phi(v^{*})(\mu_{p} -\mu_{v^{*}} \frac{\sigma_{v^{*},p}}{\sigma_{v^{*}}^2}) dv^{*} + \frac{a}{P[V^{\ast} < a]} \int_{-\infty}^{a} \phi(v^{*}) v^{*} \frac{\sigma_{v^{*},p}}{\sigma_{v^{*}}^2} dv^{*} \end{equation}
Next, we introduce: \begin{equation}\alpha= (a-\mu_{v^{*}})/\sigma_{v}\end{equation}. We know that the integral of a pdf is its cdf, so:
\begin{equation} E[VP|V^{*}<a]=\frac{a}{\Phi(\alpha)}(\mu_{p} -\mu_{v^{*}} \frac{\sigma_{v^{*},p}}{\sigma_{v^{*}}^2}) \Phi(\alpha) + \frac{a}{\Phi(\alpha)} \int_{-\infty}^{a} \phi(v^{*}) v^{*} \frac{\sigma_{v^{*},p}}{\sigma_{v^{*}}^2} dv^{*} \end{equation}
\begin{equation} E[VP|v^{*}<a]=a(\mu_{p} -\mu_{v^{*}} \frac{\sigma_{v^{*},p}}{\sigma_{v^{*}}^2}) + \frac{a}{\Phi(\alpha)} \frac{\sigma_{v^{*},p}}{\sigma_{v^{*}}^2} \int_{-\infty}^{a} \phi(v^{*}) v^{*} dv^{*} \end{equation}
The right integral is nothing but the univariate conditional truncated expectation.
\begin{equation} E[VP|V^{*}<a]=a(\mu_{p} -\mu_{v^{*}} \frac{\sigma_{v^{*},p}}{\sigma_{v^{*}}^2}) + \frac{a}{\Phi(\alpha)} \frac{\sigma_{v^{*},p}}{\sigma_{v^{*}}^2} \Phi(\alpha)[\mu_{v^{*}}-\sigma_{v^{*}}\Phi(\alpha) \frac{\phi(\alpha)}{\Phi(\alpha)}] \end{equation}
Expanded: \begin{equation} E[VP|V^{*}<a]=a [\mu_{p} -\mu_{v^{*}} \frac{\sigma_{v^{*},p}}{\sigma_{v^{*}}^2} + \frac{\sigma_{v^{*},p}}{\sigma_{v^{*}}^2}\mu_{v^{*}} -\frac{\sigma_{v^{*},p}\phi(\alpha)}{\Phi(\alpha)\sigma^2}] \end{equation}
This looks like an elegant solution (middle terms cancel out), however when simulating I get different outcomes numerically when mean and variance of v* are not 0 and 1 respectively. Meaning I still have made an error somewhere along the way. Please note that the variance of P apparently doesn't matter (it appears not to matter when I simulate,), which seems counterintuititive but ok.
[thx Daniel for pointing out an earlier error, which I made twice in my derivation].
See code for a quick comparison:
#simulate bivariate normal
import numpy as np
import pandas as pd
from scipy.stats import norm
nsim=500000
m1=1
m2=0 #this is censored
s1= 1
s2=1
s12=0.5
a=-.3
mu=[m1, m2]
sigma=[[s1, s12],[s12, s2]]
simulation = pd.DataFrame(np.random.multivariate_normal(mu,sigma,nsim))
simulation.columns=['X','V*']
alpha= (a-m2)/s2
simulation['V']=simulation['V*']
simulation['censor']=simulation['V*']<a
simulation.loc[simulation['censor'] , 'V'] = a
#plt.scatter(simulatie['X'], simulatie['V'])
PhiA=norm.cdf(alpha, loc=m2, scale=s2)
phiA=norm.pdf(alpha, loc=m2, scale=s2)
simulation['VX']=simulation['V'].multiply(simulation['X'])
analytic_expectation=a* (m1- (s12*phiA)/(PhiA*s2*s2))
numerical_average=(simulation['VX'].multiply(simulation['censor'])).mean() *(nsim/simulation['censor'].sum())
print("numerical_average: ", numerical_average )
print("analytic_expectation: ", analytic_expectation )