diff --git a/src/NowcastAutoGP.jl b/src/NowcastAutoGP.jl
index f4072ad..e3ba8bf 100644
--- a/src/NowcastAutoGP.jl
+++ b/src/NowcastAutoGP.jl
@@ -8,7 +8,7 @@ using ProgressMeter: @showprogress, Progress, next!
 
 const GPModel = AutoGP.GPModel # re-exporting for convenience
 const GPConfig = AutoGP.GP.GPConfig # re-exporting for convenience
-export TData, GPModel, GPConfig, RandomWalk
+export TData, GPModel, GPConfig, RandomWalk, IntegratedBrownianMotion
 export create_transformed_data, get_transformations, make_and_fit_model, forecast,
     forecast_with_nowcasts, create_nowcast_data
 
diff --git a/src/extension_kernels.jl b/src/extension_kernels.jl
index e6c523d..cccbb35 100644
--- a/src/extension_kernels.jl
+++ b/src/extension_kernels.jl
@@ -53,3 +53,75 @@ function rescale(node::RandomWalk, t::LinearTransform)
 end
 
 pretty(node::RandomWalk) = @sprintf("RW(%1.2f; %1.2f)", node.origin, node.amplitude)
+
+##########################################
+### IntegratedBrownianMotion kernel ######
+##########################################
+
+@doc raw"""
+    IntegratedBrownianMotion(origin[, amplitude=1])
+
+Once-integrated Brownian motion (integrated Wiener process) covariance kernel.
+
+```math
+k(t, t') = \theta_2 \, \frac{a^2 (3b - a)}{6},
+\qquad a = \min(t, t') - \theta_1, \quad b = \max(t, t') - \theta_1
+```
+
+This is the covariance of ``X(t) = \int_{\theta_1}^{t} W(s)\,\mathrm{d}s``, the integral of a
+Brownian motion that starts at the `origin` ``\theta_1`` (where both the value and the
+variance are zero); `amplitude` ``\theta_2`` scales the variance. Draws are smoother
+(once-differentiable) than [`RandomWalk`](@ref) draws — an integrated-random-walk prior, a
+natural choice for trends whose *rate of change* drifts like a random walk.
+
+The kernel is positive-definite only when `origin <= min(t)` over the evaluated time points,
+so that every `a, b >= 0`.
+
+Like [`RandomWalk`](@ref), this kernel is defined in `NowcastAutoGP` (not `AutoGP`) by
+extending `AutoGP`'s `GP` interface (`eval_cov`, `reparameterize`, `rescale`).
+
+### Connection to cubic splines
+
+The predictive mean of a GP with kernel composition of primitives `Const` + `Linear` + `IntegratedBrownianMotion`
+is a cubic spline (piecewise cubic polynomial) between data points and linear outside the data range, see [Rasmussen & Williams, 2006, §6.3].
+This means that a GP with this kernel composition represents a Bayesian cubic spline model, with the `IntegratedBrownianMotion` kernel controlling the smoothness of the spline.
+Note that the kernel covariance here seems to differ in form from the one in Rasmussen & Williams, but noting that |t - t'| = b - a, the two forms are equivalent.
+"""
+struct IntegratedBrownianMotion <: LeafNode
+    origin::Real
+    amplitude::Real
+    IntegratedBrownianMotion(origin::Real, amplitude::Real = 1) = new(origin, amplitude)
+end
+
+function eval_cov(node::IntegratedBrownianMotion, t1, t2)
+    a = min(t1, t2) - node.origin
+    b = max(t1, t2) - node.origin
+    return node.amplitude * a^2 * (3b - a) / 6
+end
+
+function eval_cov(node::IntegratedBrownianMotion, ts::Vector{Float64})
+    a = min.(ts, ts') .- node.origin
+    b = max.(ts, ts') .- node.origin
+    return node.amplitude .* a .^ 2 .* (3 .* b .- a) ./ 6
+end
+
+# Input transform f(t) = a*t + b (a = t.slope > 0, b = t.intercept). The covariance is a
+# degree-3 homogeneous function of (t - origin, t' - origin) (a double integral of min, which
+# is degree-1), so under f the origin transforms as for `RandomWalk` while the amplitude
+# absorbs a factor of slope^3:
+#   origin'    = (origin - b) / a
+#   amplitude' = a^3 * amplitude
+function reparameterize(node::IntegratedBrownianMotion, t::LinearTransform)
+    origin = (node.origin - t.intercept) / t.slope
+    amplitude = t.slope^3 * node.amplitude
+    return IntegratedBrownianMotion(origin, amplitude)
+end
+
+# Output transform Y = a*X + b scales variance by a^2; the origin is unchanged.
+function rescale(node::IntegratedBrownianMotion, t::LinearTransform)
+    return IntegratedBrownianMotion(node.origin, t.slope^2 * node.amplitude)
+end
+
+function pretty(node::IntegratedBrownianMotion)
+    return @sprintf("IBM(%1.2f; %1.2f)", node.origin, node.amplitude)
+end
diff --git a/test/test_integrated_brownian_motion_kernel.jl b/test/test_integrated_brownian_motion_kernel.jl
new file mode 100644
index 0000000..4bf0e53
--- /dev/null
+++ b/test/test_integrated_brownian_motion_kernel.jl
@@ -0,0 +1,83 @@
+@testitem "IntegratedBrownianMotion eval_cov (scalar form)" begin
+    using AutoGP
+
+    node = IntegratedBrownianMotion(0.0, 2.0)
+    # For origin = 0 and s <= t the covariance is amplitude * (s^2 t / 2 - s^3 / 6).
+    s, t = 0.3, 0.5
+    @test AutoGP.GP.eval_cov(node, s, t) ≈ 2.0 * (s^2 * t / 2 - s^3 / 6)
+    @test AutoGP.GP.eval_cov(node, t, s) ≈ AutoGP.GP.eval_cov(node, s, t)   # symmetric
+    @test AutoGP.GP.eval_cov(node, 0.0, 0.4) == 0.0                         # zero variance at origin
+
+    # a non-zero origin shifts time: only (t - origin) enters the kernel
+    shifted = IntegratedBrownianMotion(-1.0, 1.0)
+    a, b = 0.3 - (-1.0), 0.5 - (-1.0)
+    @test AutoGP.GP.eval_cov(shifted, 0.3, 0.5) ≈ a^2 * (3b - a) / 6
+
+    # default amplitude is 1
+    @test IntegratedBrownianMotion(0.0).amplitude == 1
+end
+
+@testitem "IntegratedBrownianMotion eval_cov (matrix form) matches scalar and is PSD" begin
+    using AutoGP
+    using Random: Random
+
+    node = IntegratedBrownianMotion(0.0, 1.5)
+    ts = collect(0.0:0.25:1.0)
+    C = AutoGP.GP.eval_cov(node, ts)
+
+    @test size(C) == (length(ts), length(ts))
+    # matrix entries equal the pairwise scalar evaluation
+    for i in eachindex(ts), j in eachindex(ts)
+        @test C[i, j] == AutoGP.GP.eval_cov(node, ts[i], ts[j])
+    end
+    @test C == C'                                   # symmetric
+
+    # positive semi-definite (origin = 0 <= min(ts)): vᵀ C v >= 0 for all v.
+    # Checked without LinearAlgebra via the quadratic form on many random vectors.
+    rng = Random.MersenneTwister(1234)
+    quad(v) = sum(v[i] * C[i, j] * v[j] for i in eachindex(v), j in eachindex(v))
+    for _ in 1:1000
+        @test quad(randn(rng, length(ts))) >= -1.0e-10
+    end
+end
+
+@testitem "IntegratedBrownianMotion reparameterize matches input warping" begin
+    using AutoGP
+
+    node = IntegratedBrownianMotion(-0.5, 2.0)
+    t = AutoGP.GP.LinearTransform(3.0, -1.0)   # f(x) = 3x - 1, slope > 0
+    warped = AutoGP.GP.reparameterize(node, t)
+
+    # defining property: k(reparameterize(n, t), s, u) == k(n, f(s), f(u))
+    for (s, u) in ((0.2, 0.7), (0.5, 0.5), (0.9, 0.1))
+        fs = t.slope * s + t.intercept
+        fu = t.slope * u + t.intercept
+        @test AutoGP.GP.eval_cov(warped, s, u) ≈ AutoGP.GP.eval_cov(node, fs, fu)
+    end
+
+    # closed-form parameters (origin like RandomWalk; amplitude is degree-3 in slope)
+    @test warped.origin ≈ (node.origin - t.intercept) / t.slope
+    @test warped.amplitude ≈ t.slope^3 * node.amplitude
+end
+
+@testitem "IntegratedBrownianMotion rescale matches output scaling" begin
+    using AutoGP
+
+    node = IntegratedBrownianMotion(-0.5, 2.0)
+    t = AutoGP.GP.LinearTransform(4.0, 7.0)   # output warp Y = 4X + 7
+    scaled = AutoGP.GP.rescale(node, t)
+
+    # output scaling multiplies the variance by slope^2; origin is unchanged
+    @test scaled.origin == node.origin
+    @test scaled.amplitude ≈ t.slope^2 * node.amplitude
+    for (s, u) in ((0.2, 0.7), (0.5, 0.5))
+        @test AutoGP.GP.eval_cov(scaled, s, u) ≈ t.slope^2 * AutoGP.GP.eval_cov(node, s, u)
+    end
+end
+
+@testitem "IntegratedBrownianMotion pretty printing" begin
+    using AutoGP
+
+    @test AutoGP.GP.pretty(IntegratedBrownianMotion(0.0, 2.0)) == "IBM(0.00; 2.00)"
+    @test AutoGP.GP.pretty(IntegratedBrownianMotion(-1.5, 0.5)) == "IBM(-1.50; 0.50)"
+end