  We present recent progress towards a global determination of parton
distribution functions (PDFs) at approximate N$^3$LO (aN$^3$LO) accuracy within
the NNPDF framework. We construct a parametrisation of the
$\mathcal{O}(\alpha_s^4)$ QCD splitting functions and anomalous dimensions
reproducing all known exact results, estimate the associated missing and
incomplete higher order uncertainties (MHOU and IHOU, respectively), and
implement it in the open-source DGLAP code EKO enabling PDFs to be evolved at
aN$^3$LO accuracy in the NNPDF fitting framework. We compare aN$^3$LO
calculation of splitting functions with the results of lower perturbative
orders and quantify the impact of the various sources of theoretical
uncertainties.
