We consider the computation by simulation and neural net regression of conditional expectations, or more general elicitable statistics, of functionals of processes $(X, Y )$. Here an exogenous component $Y$ (Markov by itself) is time-consuming to simulate, while the endogenous component $X$ (jointly Markov with $Y$) is quick to simulate given $Y$, but is responsible for most of the variance of the simulated payoff. To address the related variance issue, we introduce a conditionally independent, hierarchical simulation scheme, where several paths of $X$ are simulated for each simulated path of $Y$. We analyze the statistical convergence of the regression learning scheme based on such block-dependent data. We derive heuristics on the number of paths of $Y$ and, for each of them, of $X$, that should be simulated. The resulting algorithm is implemented on a graphics processing unit (GPU) combining Python/CUDA and learning with PyTorch. A CVA case study with a nested Monte Carlo benchmark shows that the hierarchical simulation technique is key to the success of the learning approach., Comment: This article has been accepted for publication in Mathematical Finance, published by Wiley