Under an appropriate regular variation condition, the affinely normalized partial sums of a sequence of independent and identically distributed random variables converges weakly to a non-Gaussian stable random variable. A functional version of this is known to be true as well, the limit process being a stable Lévy process. The main result in the paper is that for a stationary, regularly varying sequence for which clusters of high-threshold excesses can be broken down into asymptotically independent blocks, the properly centered partial sum process still converges to a stable Lévy process. Due to clustering, the Lévy triple of the limit process can be different from the one in the independent case. The convergence takes place in the space of càdlàg functions endowed with Skorohod's $M_1$ topology, the more usual $J_1$ topology being inappropriate as the partial sum processes may exhibit rapid successions of jumps within temporal clusters of large values, collapsing in the limit to a single jump. The result rests on a new limit theorem for point processes which is of independent interest. The theory is applied to moving average processes, squared GARCH(1,1) processes and stochastic volatility models.