Alright, with new design of DynamicalSystems.jl v3 it is rather straightforward to compute the maximum lyapunov exponent of arbitrary dynamical systems. Below, I am attaching a code that computes the Lyapunov exponent of the duffing oscillator first as a normal system and then as a stroboscopic map:

using ChaosTools

# %% Stroboscopic
function duffing_rule(x, p, t)
    ω, f, d, β = p
    dx1 = x[2]
    dx2 = f*cos(ω*t) - β*x[1] - x[1]^3 - d * x[2]
    return SVector(dx1, dx2)
end

u0 = [0.1, 0.25]
p0 = [1.0, 0.3, 0.2, -1]
T = 2π/1.0

duffing_raw = CoupledODEs(duffing_rule, u0, p0)
duffing_map = StroboscopicMap(duffing_raw, T)

λspec = lyapunovspectrum(duffing_raw, 10000)
λmax = lyapunov(duffing_raw, 10000)
λmax_smap = lyapunov(duffing_map, 10000)

@show λspec
@show λmax, λmax*T, λmax_smap

The output is:

λspec = [0.15894265996257015, -0.35894237397571577]
(λmax, λmax * T, λmax_smap) = (0.15884881282285493, 0.9980765267914823, 2.3432028585153932)

Now, as you see above, the exponent we get from the stroboscopic map does not match the period times the exponent of the normal system. Theory says it "should". I am attaching here two pages form the book of Politi and Pikovsky

Lyapunov_exponents.pdf

where they say that for a stroboscopic map one expects that the Lyapunov exponent is the period times the exponent of the normal time system.

However, I've checked the code extensively and I am not sure there is a mistake. The code does what its supposed to do, treating the stroboscopic map as a deterministic iterated map. I hope there is some easy to find bug somewhere so that we don't have to dig deep into theory to check if the statement of the book is actually correct..