add RW primitive kernel

Project.toml

@@ -9,6 +9,7 @@ BoxCox = "1248164d-f7a6-4bdb-8e8d-8c4a187b3ce6"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 
 [compat]
@@ -17,5 +18,6 @@ BoxCox = "0.3.7"
 Dates = "1.10"
 LinearAlgebra = "1.10"
 LogExpFunctions = "0.3.29, 1"
+Printf = "1.10"
 ProgressMeter = "1.10"
 julia = "1.11"

src/NowcastAutoGP.jl

@@ -3,16 +3,18 @@ using AutoGP, Dates
 using BoxCox: BoxCoxTransformation, confint, fit
 using LinearAlgebra: BLAS
 using LogExpFunctions: logit, logistic
+using Printf: @sprintf
 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
+export TData, GPModel, GPConfig, RandomWalk
 export create_transformed_data, get_transformations, make_and_fit_model, forecast,
     forecast_with_nowcasts, create_nowcast_data
 
 include("transformations.jl")
 include("TData.jl")
+include("extension_kernels.jl")
 include("make_and_fit_model.jl")
 include("create_nowcast_data.jl")
 include("forecasting.jl")

src/extension_kernels.jl

@@ -0,0 +1,55 @@
+import AutoGP.GP: LeafNode, eval_cov, reparameterize, rescale, pretty, LinearTransform
+
+#############################
+### RandomWalk kernel #######
+#############################
+
+@doc raw"""
+    RandomWalk(origin[, amplitude=1])
+
+Random walk (Wiener process) covariance kernel.
+
+```math
+k(t, t') = \theta_2 \left( \min(t, t') - \theta_1 \right)
+```
+
+The `origin` ``\theta_1`` is the time at which the process starts (where the variance is
+zero); `amplitude` ``\theta_2`` scales the per-unit-time variance. Draws from this kernel are
+continuous-time random walks (Brownian motion): non-stationary, with variance growing linearly
+away from `origin`.
+
+The kernel is positive-definite only when `origin <= min(t)` over the evaluated time points,
+so that every `min(t, t') - origin >= 0`. This is the GP analogue of an i.i.d.-increment
+random walk.
+
+This kernel is defined in `NowcastAutoGP` (not `AutoGP`) by extending `AutoGP`'s `GP` interface
+(`eval_cov`, `reparameterize`, `rescale`). It mirrors AutoGP's primitive kernel structure.
+"""
+struct RandomWalk <: LeafNode
+    origin::Real
+    amplitude::Real
+    RandomWalk(origin::Real, amplitude::Real = 1) = new(origin, amplitude)
+end
+
+eval_cov(node::RandomWalk, t1, t2) = node.amplitude * (min(t1, t2) - node.origin)
+
+function eval_cov(node::RandomWalk, ts::Vector{Float64})
+    return node.amplitude .* (min.(ts, ts') .- node.origin)
+end
+
+# Input transform f(t) = a*t + b (a = t.slope > 0, b = t.intercept). We need
+# k(f(t), f(u); θ) = k(t, u; θ'). Since min(a*t+b, a*u+b) = a*min(t,u) + b,
+#   origin'    = (origin - b) / a       (identical to Linear's intercept rule)
+#   amplitude' = a * amplitude          (degree-1 in t, so scales linearly with slope)
+function reparameterize(node::RandomWalk, t::LinearTransform)
+    origin = (node.origin - t.intercept) / t.slope
+    amplitude = t.slope * node.amplitude
+    return RandomWalk(origin, amplitude)
+end
+
+# Output transform Y = a*X + b scales variance by a^2; the origin is unchanged.
+function rescale(node::RandomWalk, t::LinearTransform)
+    return RandomWalk(node.origin, t.slope^2 * node.amplitude)
+end
+
+pretty(node::RandomWalk) = @sprintf("RW(%1.2f; %1.2f)", node.origin, node.amplitude)

test/test_random_walk_kernel.jl

@@ -0,0 +1,81 @@
+@testitem "RandomWalk eval_cov (scalar form)" begin
+    using AutoGP
+
+    node = RandomWalk(0.0, 2.0)
+    # k(t, t') = amplitude * (min(t, t') - origin)
+    @test AutoGP.GP.eval_cov(node, 0.3, 0.5) == 2.0 * 0.3
+    @test AutoGP.GP.eval_cov(node, 0.5, 0.3) == 2.0 * 0.3   # symmetric
+    @test AutoGP.GP.eval_cov(node, 0.4, 0.4) == 2.0 * 0.4   # diagonal
+
+    # a non-zero origin shifts the covariance by amplitude * origin
+    shifted = RandomWalk(-1.0, 1.0)
+    @test AutoGP.GP.eval_cov(shifted, 0.3, 0.5) == (0.3 - (-1.0))
+
+    # default amplitude is 1
+    @test RandomWalk(0.0).amplitude == 1
+end
+
+@testitem "RandomWalk eval_cov (matrix form) matches scalar and is PSD" begin
+    using AutoGP
+    using Random: Random
+
+    node = RandomWalk(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 "RandomWalk reparameterize matches input warping" begin
+    using AutoGP
+
+    node = RandomWalk(-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 Linear; amplitude is degree-1 in slope)
+    @test warped.origin ≈ (node.origin - t.intercept) / t.slope
+    @test warped.amplitude ≈ t.slope * node.amplitude
+end
+
+@testitem "RandomWalk rescale matches output scaling" begin
+    using AutoGP
+
+    node = RandomWalk(-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 "RandomWalk pretty printing" begin
+    using AutoGP
+
+    @test AutoGP.GP.pretty(RandomWalk(0.0, 2.0)) == "RW(0.00; 2.00)"
+    @test AutoGP.GP.pretty(RandomWalk(-1.5, 0.5)) == "RW(-1.50; 0.50)"
+end