This paper introduces a novel simulation smoothing method for state space models. The method can be used to compute smoothed estimates of the states and nonlinear functions of the states, and it allows for visualization of the joint smoothing distribution. The simulation smoother is based on the extremum Monte Carlo method. It uses simulated data from the model to estimate the conditional density functions from the backward decomposition of the joint smoothing density. The approach is generally applicable and deals naturally with missing data, as well as measurements that become available at mixed frequencies. The method is illustrated via examples with missing data, multimodal distributions, and intractable model densities. The flexibility and computational efficiency of the approach is demonstrated in an empirical application to a time series of Bitcoin based on the stochastic volatility model with stable errors.