I am trying to understand how fitting a model using MCMC works. Is there a loss function that is optimized?

Or is it simply a case of more draws from the distribution amount to a more complete description of the posterior and therefore more accurate parameters?

In particular I am referring to how a BSTS time series model is fitted using MCMC draws as described in Scott and Varian "Predicting the Present with Bayesian Structural Time Series", 2013.

In that paper a state space model of a time series is described:

$y_t = Z_t^T \alpha_t + \epsilon_t $.

$\alpha_{t+1} = T_t \alpha_t + R_t\eta_t$

Let

$\theta = (Z, T, R, \epsilon, \eta)$

and

$\textbf{y} = y_1,....y_n$ their time series data.

They then use MCMC to simulate the posterior distribution of the parameters of the model given their data $p(\theta|\textbf{y})$ - and presumably (they don't actually states this, I am assuming) they select $\theta$ that maximizes $p(\theta|\textbf{y})$.

In R the BSTS model takes a number of MCMC draws as an input parameter - and I am trying to figure how to choose that number.

If it is the second case (more draws = better simulation of the distribution), how do we decided the number of iterations and how do we avoid overfitting?