diff --git a/Project.toml b/Project.toml
index 50607efaad..4e4186023b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -46,7 +46,7 @@ LinearAlgebra = "<0.0.1, 1"
 OffsetArrays = "1"
 PDMats = "0.10, 0.11"
 Printf = "<0.0.1, 1"
-QuadGK = "2"
+QuadGK = "2.7.0"
 Random = "<0.0.1, 1"
 SparseArrays = "<0.0.1, 1"
 SpecialFunctions = "1.2, 2"
diff --git a/src/Distributions.jl b/src/Distributions.jl
index 27d18d80ab..3f55ef9eeb 100644
--- a/src/Distributions.jl
+++ b/src/Distributions.jl
@@ -3,7 +3,7 @@ module Distributions
 using StatsBase, PDMats, StatsFuns, Statistics
 using StatsFuns: logtwo, invsqrt2, invsqrt2π
 
-import QuadGK: quadgk
+import QuadGK: quadgk, quadgk!
 import Base: size, length, convert, show, getindex, rand, vec, inv
 import Base: sum, maximum, minimum, extrema, +, -, *, ==
 import Base.Math: @horner
diff --git a/src/multivariate/jointorderstatistics.jl b/src/multivariate/jointorderstatistics.jl
index 1fbed0d1b6..04ff52efec 100644
--- a/src/multivariate/jointorderstatistics.jl
+++ b/src/multivariate/jointorderstatistics.jl
@@ -166,3 +166,196 @@ function _rand!(rng::AbstractRNG, d::JointOrderStatistics, x::AbstractVector{<:R
     end
     return x
 end
+
+# Moments
+
+## Fallbacks
+
+mean(d::JointOrderStatistics{<:ContinuousUnivariateDistribution}) = _moment(d, 1)
+function var(d::JointOrderStatistics{<:ContinuousUnivariateDistribution})
+    return _moment(d, 2, mean(d))
+end
+
+function _moment(
+    d::JointOrderStatistics{<:ContinuousUnivariateDistribution},
+    p::Int,
+    μ::Union{Real,AbstractVector{<:Real}}=0;
+    rtol=sqrt(eps(promote_type(partype(d), eltype(μ)))),
+)
+    # assumes if p == 1, then μ=0 and if p>1, then μ==mean(d)
+    T = float(typeof(one(Base.promote_type(partype(d), eltype(μ)))))
+    n = d.n
+    ranks = d.ranks
+
+    μdist = mean(d.dist)
+    mp = similar(d.ranks, T)
+    logC = @. -logbeta(ranks, T(n - ranks + 1))
+    function integrand!(y, x)
+        if x < μdist
+            # for some distributions (e.g. Normal) this improves numerical stability
+            logF = logcdf(d.dist, x)
+            log1mF = log1mexp(logF)
+        else
+            log1mF = logccdf(d.dist, x)
+            logF = log1mexp(log1mF)
+        end
+        logf = logpdf(d.dist, x)
+        @. y = (x - μ)^p * exp(logC + logf + (ranks - 1) * logF + (n - ranks) * log1mF)
+        return y
+    end
+    α = eps(T) / 2
+    lower = quantile(OrderStatistic(d.dist, n, first(ranks)), α)
+    upper = quantile(OrderStatistic(d.dist, n, last(ranks)), 1 - α)
+    quadgk!(integrand!, mp, lower, upper; rtol=rtol)
+    return mp
+end
+
+## Uniform
+
+function mean(d::JointOrderStatistics{<:Uniform})
+    return d.ranks .* scale(d.dist) ./ (d.n + 1) .+ minimum(d.dist)
+end
+function var(d::JointOrderStatistics{<:Uniform})
+    n = d.n
+    ranks = d.ranks
+    return @. (scale(d.dist)^2 * ranks * (n + 1 - ranks)) / (n + 1)^2 / (n + 2)
+end
+function cov(d::JointOrderStatistics{<:Uniform})
+    n = d.n
+    ranks = d.ranks
+    c = (scale(d.dist) / (n + 1))^2 / (n + 2)
+    return broadcast(ranks, ranks') do rᵢ, rⱼ
+        rmin, rmax = minmax(rᵢ, rⱼ)
+        return rmin * (n + 1 - rmax) * c
+    end
+end
+
+## Exponential
+
+function mean(d::JointOrderStatistics{<:Exponential})
+    # Arnold, eq 4.6.6
+    θ = scale(d.dist)
+    T = float(typeof(one(θ)))
+    m = similar(d.ranks, T)
+    ks = d.n .- d.ranks
+    _harmonicnum!(m, ks)
+    Hn = _harmonicnum_from(first(m), first(ks), d.n)
+    @. m = θ * (Hn - m)
+    return m
+end
+function var(d::JointOrderStatistics{<:Exponential})
+    # Arnold, eq 4.6.7
+    θ = scale(d.dist)
+    T = float(typeof(oneunit(θ)^2))
+    v = similar(d.ranks, T)
+    ks = (d.n + 1) .- d.ranks
+    _polygamma!(v, 1, ks)
+    ϕn = _polygamma_from(1, first(v), first(ks), d.n + 1)
+    @. v = θ^2 * (v - ϕn)
+    return v
+end
+function cov(d::JointOrderStatistics{<:Exponential})
+    # Arnold, eq 4.6.8
+    v = var(d)
+    S = broadcast(d.ranks, d.ranks', v, v') do rᵢ, rⱼ, vᵢ, vⱼ
+        rᵢ < rⱼ ? vᵢ : vⱼ
+    end
+    return S
+end
+
+## Logistic
+
+function mean(d::JointOrderStatistics{<:Logistic})
+    # Arnold, eq 4.8.6
+    T = typeof(oneunit(partype(d.dist)))
+    m = H1 = similar(d.ranks, T)
+    _harmonicnum!(H1, d.n .- d.ranks)
+    if d.ranks == 1:(d.n)
+        H2 = view(H1, reverse(eachindex(H1)))
+    else
+        H2 = similar(H1)
+        _harmonicnum!(H2, d.ranks .- 1)
+    end
+    m .= scale(d.dist) .* (H2 .- H1) .+ mean(d.dist)
+    return m
+end
+function var(d::JointOrderStatistics{<:Logistic})
+    # Arnold, eq 4.8.7
+    T = typeof(oneunit(partype(d.dist))^2)
+    v = ϕ1 = similar(d.ranks, T)
+    _polygamma!(ϕ1, 1, d.ranks)
+    if d.ranks == 1:(d.n)
+        ϕ2 = view(ϕ1, reverse(eachindex(ϕ1)))
+    else
+        ϕ2 = similar(ϕ1)
+        _polygamma!(H2, 1, d.n + 1 - d.ranks)
+    end
+    v .= scale(d.dist)^2 .* (ϕ1 .+ ϕ2)
+    return v
+end
+
+## Normal
+
+function mean(d::JointOrderStatistics{<:Normal})
+    n = d.n
+    n > 5 && return _moment(d, 1)
+    return mean.(OrderStatistic.(d, d.n, d.ranks; check_args=false))
+end
+
+## AffineDistribution
+
+function mean(d::JointOrderStatistics{<:AffineDistribution})
+    σ = scale(d.dist)
+    r = σ ≥ 0 ? d.ranks : d.n + 1 .- d.ranks
+    dρ = JointOrderStatistics(d.dist.ρ, d.n, r; check_args=false)
+    return mean(dρ) .* σ .+ d.dist.μ
+end
+function var(d::JointOrderStatistics{<:AffineDistribution})
+    σ = scale(d.dist)
+    r = σ ≥ 0 ? d.ranks : d.n + 1 .- d.ranks
+    dρ = JointOrderStatistics(d.dist.ρ, d.n, r; check_args=false)
+    return var(dρ) * σ^2
+end
+function cov(d::JointOrderStatistics{<:AffineDistribution})
+    σ = scale(d.dist)
+    r = σ ≥ 0 ? d.ranks : d.n + 1 .- d.ranks
+    dρ = JointOrderStatistics(d.dist.ρ, d.n, r; check_args=false)
+    return cov(dρ) * σ^2
+end
+
+## Common utilities
+
+# assume ns are sorted in increasing or decreasing order
+function _harmonicnum!(Hns, ns::AbstractVector{<:Int})
+    Hk = zero(eltype(Hns))
+    k = 0
+    iter = if last(ns) ≥ first(ns)
+        zip(eachindex(Hns), ns)
+    else
+        Iterators.reverse(zip(eachindex(Hns), ns))
+    end
+    for (i, n) in iter
+        Hns[i] = Hk = _harmonicnum_from(Hk, k, n)
+        k = n
+    end
+    return Hns
+end
+
+# assume ns are sorted in increasing or decreasing order
+function _polygamma!(ϕns, m::Int, ns::AbstractVector{<:Int})
+    if last(ns) ≥ first(ns)
+        i = lastindex(ϕns)
+        k = last(ns)
+        iter = Iterators.reverse(zip(eachindex(ϕns), ns))
+    else
+        i = firstindex(ϕns)
+        k = first(ns)
+        iter = zip(eachindex(ϕns), ns)
+    end
+    ϕns[i] = ϕk = polygamma(m, k)
+    for (i, n) in Iterators.drop(iter, 1)
+        ϕns[i] = ϕk = _polygamma_from(m, ϕk, k, n)
+        k = n
+    end
+    return ϕns
+end
diff --git a/src/univariate/orderstatistic.jl b/src/univariate/orderstatistic.jl
index 1a7055ef91..db73a5738a 100644
--- a/src/univariate/orderstatistic.jl
+++ b/src/univariate/orderstatistic.jl
@@ -106,3 +106,266 @@ function rand(rng::AbstractRNG, d::OrderStatistic)
     b = _uniform_orderstatistic(d)
     return T(quantile(d.dist, rand(rng, b)))
 end
+
+# Moments
+
+## Fallbacks
+
+mean(d::OrderStatistic) = _moment(d, 1)
+var(d::OrderStatistic) = _moment(d, 2, _moment(d, 1))
+function skewness(d::OrderStatistic)
+    m = mean(d)
+    return _moment(d, 3, m) / _moment(d, 2, m)^(3//2)
+end
+function kurtosis(d::OrderStatistic)
+    m = mean(d)
+    return _moment(d, 4, m) / _moment(d, 2, m)^2 - 3
+end
+
+function _moment(
+    d::OrderStatistic{<:ContinuousUnivariateDistribution},
+    p::Int,
+    μ::Real=0;
+    rtol=sqrt(eps(promote_type(partype(d), typeof(μ)))),
+)
+    # assumes if p == 1, then μ=0 and if p>1, then μ==mean(d)
+    T = float(typeof(one(Base.promote_type(partype(d), typeof(μ)))))
+    n = d.n
+    rank = d.rank
+
+    μdist = mean(d.dist)
+    isfinite(μ) && isfinite(μdist) || return T(NaN)
+    if isodd(p) && isodd(n) && rank == (n + 1) ÷ 2 && _issymmetric(d.dist)
+        # for symmetric distributions, distribution of median is also symmetric, so all of
+        # its odd central moments are 0.
+        return p == 1 ? μdist : zero(T)
+    end
+
+    logC = -logbeta(rank, T(n - rank + 1))
+    function integrand(x)
+        if x < μ
+            # for some distributions (e.g. Normal) this improves numerical stability
+            logF = logcdf(d.dist, x)
+            log1mF = log1mexp(logF)
+        else
+            log1mF = logccdf(d.dist, x)
+            logF = log1mexp(log1mF)
+        end
+        return (x - μ)^p *
+               exp(logC + logpdf(d.dist, x) + (rank - 1) * logF + (n - rank) * log1mF)
+    end
+    α = eps(T) / 2
+    lower = quantile(d, α)
+    upper = quantile(d, 1 - α)
+    return first(quadgk(integrand, lower, upper; rtol=rtol))
+end
+
+function _moment(d::OrderStatistic{<:DiscreteUnivariateDistribution}, p::Int, μ::Real=0)
+    T = float(typeof(one(Base.promote_type(partype(d), typeof(μ)))))
+
+    if isbounded(d.dist)
+        xs = support(d.dist)
+    elseif eltype(d.dist) <: Integer
+        # Note: this approximation can fail in the unlikely case of an atom with huge
+        # magnitude just outside of the bulk.
+        α = eps(T) / 2
+        xmin = quantile(d, α)
+        xmax = quantile(d, 1 - α)
+        xs = xmin:xmax
+    else
+        throw(
+            ArgumentError(
+                "Moments can only be computed for bounded distributions or those with integer support.",
+            ),
+        )
+    end
+    length(xs) == 1 && p == 1 && return first(xs) - μ
+
+    b = _uniform_orderstatistic(d)
+    cx = cdf(d.dist, first(xs) - 1)
+    c = cdf(b, cx)
+    mp = zero(first(xs)) * c
+    for x in xs
+        cx += pdf(d.dist, x)
+        c_new = cdf(b, cx)
+        mp += (x - μ)^p * (c_new - c)
+        c = c_new
+    end
+
+    return mp
+end
+
+_issymmetric(d) = false
+_issymmetric(::Normal) = true
+_issymmetric(::NormalCanon) = true
+_issymmetric(::Laplace) = true
+_issymmetric(::Arcsine) = true
+_issymmetric(::TDist) = true
+_issymmetric(d::Beta) = d.α == d.β
+_issymmetric(::Biweight) = true
+_issymmetric(::Triweight) = true
+_issymmetric(::SymTriangularDist) = true
+
+## Uniform
+
+mean(d::OrderStatistic{<:Uniform}) = d.rank * scale(d.dist) / (d.n + 1) + minimum(d)
+std(d::OrderStatistic{<:Uniform}) = std(_uniform_orderstatistic(d)) * scale(d.dist)
+var(d::OrderStatistic{<:Uniform}) = var(_uniform_orderstatistic(d)) * scale(d.dist)^2
+skewness(d::OrderStatistic{<:Uniform}) = skewness(_uniform_orderstatistic(d))
+kurtosis(d::OrderStatistic{<:Uniform}) = kurtosis(_uniform_orderstatistic(d))
+
+## Exponential
+
+function mean(d::OrderStatistic{<:Exponential})
+    # Arnold, eq 4.6.6
+    θ = scale(d.dist)
+    T = float(typeof(one(θ)))
+    return θ * _harmonicdiff(T, d.n - d.rank, d.n)
+end
+function var(d::OrderStatistic{<:Exponential})
+    # Arnold, eq 4.6.7
+    θ = scale(d.dist)
+    T = float(typeof(one(θ)))
+    return θ^2 * _polygamma_diff(T, 1, d.n + 1 - d.rank, d.n + 1)
+end
+function skewness(d::OrderStatistic{<:Exponential})
+    θ = scale(d.dist)
+    T = float(typeof(one(θ)))
+    return -_polygamma_diff(T, 2, d.n + 1 - d.rank, d.n + 1) /
+           _polygamma_diff(T, 1, d.n + 1 - d.rank, d.n + 1)^(3//2)
+end
+function kurtosis(d::OrderStatistic{<:Exponential})
+    θ = scale(d.dist)
+    T = float(typeof(one(θ)))
+    return _polygamma_diff(T, 3, d.n + 1 - d.rank, d.n + 1) /
+           _polygamma_diff(T, 1, d.n + 1 - d.rank, d.n + 1)^2
+end
+
+## Logistic
+
+function mean(d::OrderStatistic{<:Logistic})
+    # Arnold, eq 4.8.6
+    T = typeof(oneunit(partype(d.dist)))
+    return scale(d.dist) * _harmonicdiff(T, d.n - d.rank, d.rank - 1) + mean(d.dist)
+end
+function var(d::OrderStatistic{<:Logistic})
+    # Arnold, eq 4.8.7
+    σ = scale(d.dist)
+    T = float(typeof(one(σ)))
+    return σ^2 * _polygamma_sum(T, 1, d.n + 1 - d.rank, d.rank)
+end
+function skewness(d::OrderStatistic{<:Logistic})
+    σ = scale(d.dist)
+    T = float(typeof(one(σ)))
+    return _polygamma_diff(T, 2, d.rank, d.n + 1 - d.rank) /
+           _polygamma_sum(T, 1, d.rank, d.n + 1 - d.rank)^(3//2)
+end
+function kurtosis(d::OrderStatistic{<:Logistic})
+    σ = scale(d.dist)
+    T = float(typeof(one(σ)))
+    return _polygamma_sum(T, 3, d.rank, d.n + 1 - d.rank) /
+           _polygamma_sum(T, 1, d.rank, d.n + 1 - d.rank)^2
+end
+
+## Normal
+
+function mean(d::OrderStatistic{<:Normal})
+    n = d.n
+    n > 5 && return _moment(d, 1)
+    rank = d.rank
+    μ = mean(d.dist)
+    σ = scale(d.dist)
+    T = float(typeof(one(σ)))
+    # Arnold §4.9
+    isodd(n) && rank == (n + 1) ÷ 2 && return μ
+    n == 2 && return (2 * rank - 3) * σ / sqrtπ + μ
+    n == 3 && return (rank - 2) * 3σ / 2 / sqrtπ + μ
+    I2 = atan(T(sqrt2)) / π
+    I3 = (6I2 - 1) / 4
+    r = max(rank, n - rank + 1)
+    s = (-1)^(d.rank != r)
+    if n == 4
+        c = 6s * σ / sqrtπ
+        r == 4 && return c * I2 + μ
+        r == 5 && return c * (1 - 3I2) + μ
+    end
+    if n == 5
+        c = 10s * σ / sqrtπ
+        r == 5 && return c * I3 + μ
+        r == 4 && return c * (3I2 - 4I3) + μ
+    end
+end
+
+## AffineDistribution
+
+function mean(d::OrderStatistic{<:AffineDistribution})
+    σ = scale(d.dist)
+    r = σ ≥ 0 ? d.rank : d.n - d.rank + 1
+    dρ = OrderStatistic(d.dist.ρ, d.n, r; check_args=false)
+    return mean(dρ) * σ + d.dist.μ
+end
+function var(d::OrderStatistic{<:AffineDistribution})
+    σ = scale(d.dist)
+    r = σ ≥ 0 ? d.rank : d.n - d.rank + 1
+    dρ = OrderStatistic(d.dist.ρ, d.n, r; check_args=false)
+    return var(dρ) * σ^2
+end
+function skewness(d::OrderStatistic{<:AffineDistribution})
+    σ = scale(d.dist)
+    r = σ ≥ 0 ? d.rank : d.n - d.rank + 1
+    return sign(σ) * skewness(OrderStatistic(d.dist.ρ, d.n, r; check_args=false))
+end
+function kurtosis(d::OrderStatistic{<:AffineDistribution})
+    r = scale(d.dist) ≥ 0 ? d.rank : d.n - d.rank + 1
+    return kurtosis(OrderStatistic(d.dist.ρ, d.n, r; check_args=false))
+end
+
+# Common utilities
+
+_harmonicnum(T::Type{<:Real}, n::Int) = _harmonicnum_from(zero(T), 0, n)
+
+function _harmonicnum_from(Hm::Real, m::Int, n::Int)
+    # m ≤ n
+    (n - m) < 25 && return sum(Base.Fix1(/, one(Hm)), (m + 1):n; init=Hm)
+    return digamma(oftype(Hm, n + 1)) + Base.MathConstants.eulergamma
+end
+
+function _harmonicdiff(T::Type{<:Real}, m::Int, n::Int)
+    # TODO: improve heuristic
+    d = n - m
+    m, n = minmax(m, n)
+    abs(d) < 50 && return sign(d) * sum(Base.Fix1(/, one(T)), (m + 1):n; init=zero(T))
+    Hm = _harmonicnum(T, m)
+    Hn = _harmonicnum_from(Hm, m, n)
+    return sign(d) * (Hn - Hm)
+end
+
+function _polygamma_from(m, ϕk::Real, k::Int, n::Int)
+    # backwards recurrence is more stable than forwards
+    gap = k - n
+    gap > 10 || gap < 0 && return polygamma(m, oftype(ϕk, n))
+    num = (-1)^(m + 1) * oftype(ϕk, factorial(m))
+    f = Base.Fix1(/, num) ∘ Base.Fix2(^, m + 1)
+    return sum(f, (k - 1):-1:n; init=ϕk)
+end
+
+function _polygamma_diff(T::Type{<:Real}, m::Int, k::Int, n::Int)
+    d = n - k
+    k, n = minmax(k, n)
+    s = -sign(d)
+    if abs(d) ≤ 10
+        num = (-1)^m * s * T(factorial(m))
+        f = Base.Fix1(/, num) ∘ Base.Fix2(^, m + 1)
+        return sum(f, k:(n - 1); init=zero(T))
+    end
+    ϕn = polygamma(m, T(n))
+    ϕk = _polygamma_from(m, ϕn, n, k)
+    return s * (ϕn - ϕk)
+end
+
+function _polygamma_sum(T::Type{<:Real}, m::Int, k::Int, n::Int)
+    k, n = minmax(k, n)
+    ϕn = polygamma(m, T(n))
+    ϕk = _polygamma_from(m, ϕn, n, k)
+    return ϕn + ϕk
+end