So my question is: is Pauli-repulsion a phenomenon that has also not
yet been explained in terms of any of the three other forces that we
know of?
$\def\ket#1{|#1\rangle} \let\up=\uparrow \let\dn=\downarrow
\def\PD#1#2{{\partial#1\over\partial#2}}$
There is no repulsion and no unexplained force. I would also add that PEP is an outdated way of describing the matter. In QM you should rather speak of antisymmetry of fermion states. It's only when we build up a many particle state as a tensor product of one particle states that antisymmetry forces us to keep only different states for each single particle. A simple example with two particles will explain this (I hope).
Two identical fermions in an infinite well
Consider particles in one dimension, constrained in a segment $0\le
x\le L$ (what is usally called an "infinite potential well"). Energy eigenfunctions (standing waves) are sinusoidal waves vanishing at boundaries:
$$\psi_n = \sin {n\,\pi\,x \over L} \qquad (n = 1,2,\ldots)$$
(these aren't normalized, but it's of no consequence for my present purposes.) The corresponding energy eigenvalues are
$$E_n = {n^2 h^2 \over 8\,m\,L^2}.\tag1$$
A short derivation follows, which you may skip with no harm.
$\psi_n$ has wavelength $2L/n$, then momentum
$$p = {h \over \lambda} = {n\,h \over 2\,L}.$$
Then energy (only kinetic) is
$$E_n = {p^2 \over 2\,m} = {n^2 h^2 \over 8\,m\,L^2}.$$
Assume your particles are non-interacting spin 1/2 fermions. Then above expression for energy eigenfunction is to be supplemented by specifying the spin state. Then Dirac's ket notation is preferable:
$$\ket{n\up} \quad \hbox{or} \quad
\ket{n\dn}$$
both belonging to $E_n$ eigenvalue.
If your system consists of just two particles, a set of base kets would be obtained by taking tensor products, which in Dirac's notation are written just putting two kets one after another. E.g.
$$\ket{m\up} \ket{n\up} \quad \ket{m\up} \ket{n\dn} \quad \ket{m\dn} \ket{n\up} \quad \ket{m\dn} \ket{n\dn}$$
for all positive integers $m$, $n$. A shorthand may be used:
$$\ket{m\up\,;\,n\up} \
\ket{m\up\,;\,n\dn} \
\ket{m\dn\,;\,n\up} \
\ket{m\dn\,;\,n\dn} \tag2$$
where labels preceding ";" refer to first particle, those following to the second.
But states in (2) are wrong for identical fermion particles, as they
aren't antisymmetrized. The right ones are
$$\eqalign{
&\ket{m\up\,;\,n\up} - \ket{n\up\,;\,m\up} \qquad
\ket{m\up\,;\,n\dn} - \ket{n\dn\,;\,m\up} \cr
&\ket{m\dn\,;\,n\up} - \ket{n\up\,;\,m\dn} \qquad
\ket{m\dn\,;\,n\dn} - \ket{n\dn\,;\,m\dn} \cr}$$
(once again I'm neglecting normalization).
Observe however that if $m=n$ first and fourth expressions are identically zero, whereas second and third are the same apart for sign, thus representing the same state. This is the mathematical form PEP assumes in QM: for $m=n$ just one state exists for two particles, for $m\ne n$ there are four.
For more particles we would proceed analogously, with a somewhat higher complication.
Let's compute pressure
First of all let me remark that not fermions alone exert a pressure when confined in a finite volume. Bosons do as well. Radiation pressure is an example, and photons are bosons. So let's compute the pressure exerted by a gas of non-interacting bosons at $0\,$K, when all particles are in the ground state (this isn't forbidden for bosons).
If we have $N$ particles, overall energy is given by (1) taken for $n=1$ and multiplied by $N$;
$$E = {N h^2 \over 8\,m\,L^2}.$$
As we are in one dimension we'll speak of force, not of pressure. It's most easily computed by
$$F = -\PD EL = {N h^2 \over 4\,m\,L^3}.\tag3$$
For those who find too abstract the above derivation I'll add a semiclassical one. In our box we have free particles bouncing back and forth between boundaries. Their momentum is $p=h/(2L)$. A particle hits one boundary (e.g. the left one) once in a time
$${2L \over v} = {2mL \over p} =
{4 m L^2 \over h}$$
and every time it exchanges with the boundary a momentum $2p$. Then the momentum exchanged per unit of time, i.e. the force, is
$$f = 2p\, {h \over 4 m L^2} =
{h^2 \over 4 m L^3}.$$
This holds for one particle. It's only left to multiply by $N$ to get (3).
Now for fermions
What's the difference? Simply that even at $0\,$K a fermion gas
doesn't have all particles in ground state. We've seen why it's
forbidden by antisymmetry. So we have the task to arrange an
antisymmetrical ket for $N$ particles, which sounds prohibitive.
Actually it's not so much so, but we'll follow a roundabout way, in
principle an approximated one but absolutely adequate to our purposes.
For each $n$ there are two states allowed, spin up and spin down. We
already saw that for $m=n=1$ and two particles only one state is
possible, wheres none is possibile for three. If we accept values 1 and
2 for $m$, $n$ we can accomodate up to four particles
$$\ket{1\up;1\dn;2\up;2\dn}$$
(to be antisymmetrized). So we see that for $N$ particles all states
from 1 to $N/2$ will be occupied, each by two particles with opposite
spins.
And now we are able to compute the energy:
$$E = 2\,\sum_{n=1}^{N/2} E_n =
2\,\sum_{n=1}^{N/2} {n^2 h^2 \over 8\,m\,L^2} =
{h^2 \over 4\,m\,L^2} \sum_{n=1}^{N/2} n^2$$
(the sum has to be multiplied by 2 since for every $n$ there are two
spin states). If $N$ is large we may approximate the sum to ${1 \over 24}\,N^3$ and get
$$E = {N^3 h^2 \over 96\,m\,L^2}.$$
As before
$$F = -\PD EL = {N^3 h^2 \over 48\,m\,L^3}.\tag4$$
You can see the difference between (3) and (4). Whereas for bosons
force is $\propto N$, for fermions it's $\propto N^3$, then much larger
if $N$ is large. Actually extremely larger for a white dwarf: try to
estimate how much is $N$ (number of electrons) for a star having Sun's
mass.
To be sure we should reason about pressure, not about force. This
requires leaving our naive 1D model for a more realistic 3D one. I'll
content myself to give the result
$$P = {(3\,\pi^2)}^{2/3}
\left(\!{\hbar^2 \over m}\!\right)\,{N \over V}^{\!5/3}.$$
The most important difference is in the dependence on $N$: $N^{5/3}$
instead of $N^3$. I can't explain its origin (it has to do with the
different accounting in 1D and in 3D for the one-particle states up to
$N/2$). I'll only say that even with the smaller exponent resulting
pressure is enough to counterbalance gravity for dwarfs of mass near
Sun's and size about Earth's.
A final comment
It should be clear that no mysterious force could account for our
results. Note that total energy of $N$ particles depends on a power of
$N$ and it would be hard to explain that with some interaction between
particles. Instead all depends on which and how many independent states
are allowed when identical particles are concerned. In a different
way for bosons against fermions and both different of the one that
would be used for classical particles.
As Feynman liked to say, this is the way things are.