  We establish sharp interior and boundary regularity estimates for solutions
to $\partial_t u - L u = f(t, x)$ in $I\times \Omega$, with $I \subset
\mathbb{R}$ and $\Omega \subset\mathbb{R}^n$. The operators $L$ we consider are
infinitessimal generators of stable L\'evy processes. These are linear nonlocal
operators with kernels that may be very singular. On the one hand, we establish
interior estimates, obtaining that $u$ is $C^{2s+\alpha}$ in $x$ and
$C^{1+\frac{\alpha}{2s}}$ in $t$, whenever $f$ is $C^{\alpha}$ in $x$ and
$C^{\frac{\alpha}{2s}}$ in $t$. In the case $f\in L^\infty$, we prove that $u$
is $C^{2s-\epsilon}$ in $x$ and $C^{1-\epsilon}$ in $t$, for any $\epsilon >
0$. On the other hand, we study the boundary regularity of solutions in
$C^{1,1}$ domains. We prove that for solutions $u$ to the Dirichlet problem the
quotient $u/d^s$ is H\"older continuous in space and time up to the boundary
$\partial\Omega$, where $d$ is the distance to $\partial\Omega$. This is new
even when $L$ is the fractional Laplacian.
