Solution (source code)

= Solution

At the fitted <hyperparameters>, define $K= k_0(t,t)$, $K_{**}=k_0(t^*,t^*)$, and $K_{*t}=k_0(t^*,t)$, with the matrices formed by pairwise kernel evaluations. Conditional independence of observation noise and future latent values gives
$$
\begin{pmatrix}y\\f_*\end{pmatrix}\sim N\left(0,\begin{pmatrix}V&K_{t*}\\K_{*t}&K_{**}\end{pmatrix}\right),\qquad V=K+R.
$$
The residual $f_*-K_{*t}V^{-1}y$ has zero <covariance> with $y$, so the joint normal law makes it independent of $y$. This derives the <conditional multivariate normal distribution> and the <Gaussian process regression posterior>:
$$
\boxed{f_*\mid y,t,t^*,\widehat{\mathcal H}\sim N(m_*,\Sigma_*),\qquad m_*=K_{*t}V^{-1}y,\qquad\Sigma_*=K_{**}-K_{*t}V^{-1}K_{t*}}.
$$
The <covariance> is a <Schur complement> and is positive semidefinite. When nonsingular, its joint density is $(2\pi)^{-M/2}|\Sigma_*|^{-1/2}\exp[-\frac12(f_*-m_*)^T\Sigma_*^{-1}(f_*-m_*)]$. Repeated phases can make it singular, in which case the displayed normal law is interpreted on its support.

This predicts the true latent <light curve>, so no future measurement <variance> is added. For noisy future observations, add their independent noise <covariance> to $\Sigma_*$. Fixing the <hyperparameters> omits their estimation uncertainty, as requested; integrating over their posterior would add it. The periodic model extrapolates by phase, even though all prediction times follow the observed times.