  We present a calculation of the pion and kaon Mellin moment $\langle x^3
\rangle$ extracted directly in lattice QCD using a three-derivative local
operator. We use one ensemble of gauge configurations with two degenerate
light, a strange and a charm quark ($N_f=2+1+1$) of maximally twisted mass
fermions with clover improvement. The ensemble reproduces a pion mass $\sim260$
MeV, and a kaon mass $\sim530$ MeV. Excited-states contamination is evaluated
using four values of the source-sink time separation within the range of
$1.12-1.67$ fm. We use an operator that is free of mixing, and apply a
multiplicative renormalization function calculated non-perturbatively. Our
results are converted to the $\overline{\rm MS}$ scheme and evolved at a scale
of 2 GeV, using three-loop expressions in perturbation theory. The final values
are $\langle x^3 \rangle_\pi^{u^+}=0.024(18)_{\rm stat}(2)_{\rm syst}$,
$\langle x^3 \rangle_K^{u^+}=0.035(6)_{\rm stat}(3)_{\rm syst}$, and $\langle
x^3 \rangle_K^{s^+}=0.075(5)_{\rm stat}(1)_{\rm syst}$, where the systematic
error is the uncertainty due to excited state contamination.
  We combine $\langle x^3 \rangle$ with the two lower moments to obtain the
ratios $\langle x^3 \rangle/\langle x \rangle$ and $\langle x^3 \rangle/\langle
x^2 \rangle$, as well as ratios between the pion and kaon moments. In addition,
we reconstruct the $x$-dependence of the pion and kaon PDFs via 2- and
3-parameter fits to our results. We find that the reconstruction is feasible
and that our lattice data favor a large $x$-dependence that falls as $(1-x)^2$
for both the pion and kaon PDFs. We integrate the reconstructed PDFs to extract
the higher moments with $4\leq n\leq 6$. Finally, we compare the pion and kaon
PDFs, as well as the ratios of their moments, to address the effect of SU(3)
flavor symmetry breaking.
