A turbidity time series, recorded every 15 minutes in the Wemyss catchment (Scotland) for more than a year, along with two covariates (stage height and rainfall), is considered. Turbidity time series have complex dynamics because they are non-linear, non-Normal, non-stationary, have a long memory, and present missing values. All these issues are tackled by applying Markov switching autoregressive models under the Bayesian paradigm. Bayesian model choice and inference are performed numerically through novel evolutionary Monte Carlo (EMC) algorithms which are better able to traverse a complex posterior surface, such that obtained for this class of models. The chosen model presents three hidden states and autoregressions of the fourth order. Through its implementation it was possible to efficiently fit the actual data, reconstruct the sequence of hidden states, restore the missing values, and classify the observations into a few regimes, providing new insight on turbidity dynamics associated with stage height and rainfall. Hence, stream time series show periods with discrete relationship between latent variables and statistical methods accounting for these improve understanding and prediction.