  This article uses a combination of three ideas from simulation to establish a
nearly optimal polynomial upper bound for the joint density of the stable
process and its associated supremum at a fixed time on the entire support of
the joint law. The representation of the concave majorant of the stable process
and the Chambers-Mallows-Stuck representation for stable laws are used to
define an approximation of the random vector of interest. An interpolation
technique using multilevel Monte Carlo is applied to accelerate the
approximation, allowing us to establish the infinite differentiability of the
joint density as well as nearly optimal polynomial upper bounds for the joint
mixed derivatives of any order.
