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

Or is it simply a case of more draws from the distribution amount to a more accurate description of the posterior and therefore better 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$ there 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?