We establish sharp interior and boundary regularity estimates for solutions to $\partial_t u - L u = f(t, x)$ in $I\times ��$, with $I \subset \mathbb{R}$ and $��\subset\mathbb{R}^n$. The operators $L$ we consider are infinitessimal generators of stable L��vy 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+��}$ in $x$ and $C^{1+\frac��{2s}}$ in $t$, whenever $f$ is $C^��$ in $x$ and $C^{\frac��{2s}}$ in $t$. In the case $f\in L^\infty$, we prove that $u$ is $C^{2s-��}$ in $x$ and $C^{1-��}$ in $t$, for any $��> 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��lder continuous in space and time up to the boundary $\partial��$, where $d$ is the distance to $\partial��$. This is new even when $L$ is the fractional Laplacian.