  We establish sharp higher-order H\"older regularity estimates up to the
boundary for solutions to equations of the form $\partial_t u-Lu=f(t,x)$ in
$I\times\Omega$ where $I\subset\mathbb{R}$, $\Omega\subset\mathbb{R}^n$ and $f$
is H\"older continuous. The nonlocal operators $L$ considered are those arising
in stochastic processes with jumps such as the fractional Laplacian
$(-\Delta)^s$, $s\in(0,1)$.
  Our main result establishes that, if $f$ is $C^\gamma$ is space and
$C^{\gamma/2s}$ in time, and $\Omega$ is a $C^{2,\gamma}$ domain, then $u/d^s$
is $C^{s+\gamma}$ up to the boundary in space and $u$ is $C^{1+\gamma/2s}$ up
the boundary in time, where $d$ is the distance to $\partial\Omega$. This is
the first higher order boundary regularity estimate for nonlocal parabolic
equations, and is new even for the fractional Laplacian in $C^\infty$ domains.
