Let $X$ and $Y$ be compact metric spaces, let $\pi_X \colon X \times Y \to X$ be the first coordinate projection, let $\mu$ be a Borel probability measure on $X$ of full support, let $\lambda$ be a Borel probability measure on $Y$ of full support, and let
$$ f \colon X \times Y \to X \times Y $$
be a homeomorphism. Suppose that for every continuous function $g \colon X \to \mathbb{R}$,
$$ \mu \otimes \lambda \left( (x,y) \in X \times Y \, : \, \frac{1}{N} \sum_{n=0}^{N-1} g(\pi_X(f^n(x,y))) \to \int_X g \, d\mu \ \text{ as } \, N \to \infty \right) = 1. $$
Does it follow that there exists an $f$-invariant probability measure on $X \times Y$ whose $X$-marginal is $\mu$?

If so, does there necessarily exist a Borel probability measure $\nu$ on $X \times Y$ with $\lambda$-absolutely continuous $Y$-projection such that $\frac{1}{N} \sum_{n=0}^{N-1} f^n_\ast\nu$ has a subsequence converging weakly to a measure whose $X$-marginal is $\mu$?
_______________
**Motivation.**

In the case that $f$ has a skew-product structure $f(x,y)=(\theta(x),\varphi_x(y))$ [in which case the measure $\lambda$ is irrelevant for the first question], an affirmative answer to the first question is given by Corollary 6.13 of [Hans Crauel, Random Probability Measures on Polish Spaces](https://books.google.co.uk/books/about/Random_Probability_Measures_on_Polish_Sp.html?id=RhT_aZGqmQEC&redir_esc=y) - although quoting that result is probably overkill, and I think an affirmative answer to both questions is provided (with $\nu=\mu \otimes \lambda$) by just slightly adapting the regular proof of the Krylov-Bogolyubov theorem applied to $f$. In terms of application, typically this skew-product setup concerns the scenario that $X$ is the state space of an unknown "noise" and $Y$ is the "visible" state space.

In what I'm asking about, the scenario is kind-of reversed: I'm thinking of $X$ as the visible state space, and $Y$ is a "invisible" space but is equipped with a "natural" equivalence class of measures [as represented by $\lambda$] which - in generalisation of a Riemannian volume measure on a smooth manifold - heuristically defines "experimentally accessible" sets of initial conditions. There is no skew-product assumption, meaning that there is allowed to be bi-directional feedback between the visible $X$-variable and the invisible $Y$-variable. The question is then: if I observe ergodic statistics for the visible variable, does this imply that the system has an underlying (asymptotically) stationary dynamics for the process as a whole?