There are various ways to choose the covariance kernel $k(\cdot, \cdot; \theta)$ and the hyperparameters $\theta$. Hence, the question per se is probably a bit too broad to be answered reasonably. However, I figured, I could give you an example of a class of covariance kernels that appears to be popular in Gaussian process regression and explain the influence of the hyper parameters for those.
Matérn class covariance kernels are given by $$k(x_1,x_2;\ell, \nu, \sigma) = \sigma^2\frac{2^{1-\nu}}{\Gamma(\nu)}\Bigg(\sqrt{2\nu}\frac{|x_1-x_2|}{\ell}\Bigg)^\nu K_\nu\Bigg(\sqrt{2\nu}\frac{|x_1-x_2|}{\ell}\Bigg),$$ where $K_\nu$ is the modified Bessel function of the second kind. Hence, the hyper parameter $\theta := (\ell, \nu, \sigma)$.
$\ell$ is the correlation length. If $\ell$ is small, the process at two points far apart is barely correlated. If $\ell$ is large, the process at two points far apart will still be (highly) correlated. Hence, if you assume that the function has a lot of variation in a small area of the domain, the correlation length should be rather short; otherwise long.
$\nu$ is the smoothness parameter. Samples from the prior will be $ceil(\nu)-1$ times differentiable. Hence, if you know that the function you try to determine is only continuous, but not differentiable, you should choose $\nu < 1$ and larger, if you assume more regularity. $\nu = 1/2$ corresponds to the exponential covariance, $\nu = \infty$ is the Gaussian covariance and produces samples that are infinitely often differentiable.
Each point of the Gaussian process is a normally distributed random variable. Under the prior, it has distribution $\mathrm{N}(0, \sigma^2)$. Hence, the standard deviation $\sigma$ controls the pointwise variance.
Concerning the joint likelihood. I am not really sure, what you mean. In my understanding the likelihood would be something like $p(Y|f)$. This is a function that relates the likelihood of the data to the parameter in the model that is supposed to be estimated. In this case, the parameter is the function $f$. If $\epsilon \sim \mathrm{N}(0, \gamma^2)$, the likelihood is $$ p(Y|f) = \exp\left(-\frac{1}{2 \epsilon^2} \sum_{i=1}^I (y_i - f(x_i))^2 \right).$$ The other parameters are hyperparameters of the prior and do not need to appear in the likelihood in this particular case.
