The Sun's motion relative to the Solar System barycenter is dominated by the four planets Jupiter, Saturn, Uranus and Neptune. Neptune is a lot lighter than Jupiter but it's also a lot farther away, and contrary to intuition the farther out it is the larger the induced motion in the Sun around the barycenter.
*To first order we can treat the Sun's response to each planet as circular, and the final result as the linear superposition of four circular motions.
I lack the math background to derive the answer analytically, though I am confident someone will soon.
Approximate distances of the Sun to the mathematical barycenter of each binary Sun + planet pair, "amplitude":
body mass (Sun=1) a (Rsun=1) amplitude (Rsun=1)
Sun 1.000 -
Jupiter 9.548E-04 1.119E+03 1.068
Saturn 2.858E-04 2.060E+03 0.589
Uranus 4.366E-05 4.133E+03 0.180
Neptune 5.152E-05 6.465E+03 0.333
So the interesting thing is that if it were only for Jupiter, the Sun-Jupiter barycenter would always be slightly outside the Sun.
Let's run a simulation of four circular displacements with the amplitudes as shown and periods (not that it matters, as long as they are not rational number relationships) that scale as $a^{3/2}$ which conveniently will convert my four digit numbers into irrational numbers for me. :-)
As expected based on Jupiter alone putting the barycenter just outside the Sun, when all four main planets are taken into account the barycenter spends about 54% of the time outside the Sun (or the Sun spends 54% of the time further than one of it's radii from the barycenter).
median r: 1.209
mean r: 1.189
max r: 2.17
fraction: 0.542
The histogram of radial distances shows some interesting structure, which I think should be pursued further both mathematically and numerically.
Since this is a statistical analysis (there are no harmonic relationships between planets in this model) time could just as well be sampled randomly rather than using equal spaces.

In the plot above the $x$ value for each individual two-body motion is shown with color lines, but the absolute radial distance $r$ is shown in black for the combined motion.
import numpy as np
import matplotlib.pyplot as plt
amplitude = np.array([1.068, 0.589, 0.180, 0.333])
period = np.array([1.119, 2.060, 4.133, 6.465])**1.5
time = np.arange(0, 100000, 0.01)
xd = amplitude[:, None] * np.cos(2 * np.pi * time / period[:, None])
yd = amplitude[:, None] * np.sin(2 * np.pi * time / period[:, None])
r = np.sqrt(xd.sum(axis=0)**2 + yd.sum(axis=0)**2)
bins = 0.01 * np.arange(int(100 * amplitude.sum())+10)
a, b = np.histogram(r, bins=bins)
r_median = np.median(r)
r_max = r.max()
r_mean = r.mean()
fraction = (r > 1).sum() / r.sum()
print('median r: ', r_median)
print('mean r: ', r_mean)
print('max r: ', r_max)
print('fraction: ', fraction)
n = 10000
lw = 0.5
fig, (ax1, ax2) = plt.subplots(2, 1)
for x in xd:
ax1.plot(time[:n], x[:n], linewidth=lw)
ax1.plot(time[:n], r[:n], '-k', linewidth=lw)
ax2.plot(b[1:], a)
ax2.set_title('fraction > 1 is ' + str(round(fraction, 3)))
plt.show()