Extend volatility scaling for SeparableHjmModels

 - Rename HHfInv to more general scaling_matrix
 - Add BenchmarkTimesScaling enumeration with values
   ForwardRateScaling, ZeroRateScaling, DiagonalScaling
 - Amend benchmark_times_scaling method

src/models/futures/MarkovFutureModels.jl

@@ -30,6 +30,7 @@ end
         sigma_f::BackwardFlatVolatility,
         correlation_holder::Union{CorrelationHolder, Nothing},
         quanto_model::Union{AssetModel, Nothing},
+        scaling_type::BenchmarkTimesScaling = ForwardRateScaling,
         )
 
 Create a Gausian Markov model for Future prices.
@@ -41,9 +42,10 @@ function markov_future_model(
     sigma_f::BackwardFlatVolatility,
     correlation_holder::Union{CorrelationHolder, Nothing},
     quanto_model::Union{AssetModel, Nothing},
+    scaling_type::BenchmarkTimesScaling = ForwardRateScaling,
     )
     #
-    hjm_model = gaussian_hjm_model(alias, delta, chi, sigma_f, correlation_holder, quanto_model)
+    hjm_model = gaussian_hjm_model(alias, delta, chi, sigma_f, correlation_holder, quanto_model, scaling_type)
     # manage aliases different to GaussianHjmModel
     state_alias  = [ alias * "_x_" * string(k) for k in 1:length(delta()) ]
     factor_alias = [ alias * "_f_" * string(k) for k in 1:length(delta()) ]  # ensure consistency with HJM correlation calculation
@@ -149,7 +151,7 @@ function Theta(
     @assert isnothing(X) == !state_dependent_Theta(m)
     y(t) = func_y(m.hjm_model, t)
     # make sure we do not apply correlations twice in quanto adjustment!
-    sigma_T_hyb(u) = m.hjm_model.sigma_T.HHfInv .* reshape(m.hjm_model.sigma_T.sigma_f(u), (1,:))
+    sigma_T_hyb(u) = m.hjm_model.sigma_T.scaling_matrix .* reshape(m.hjm_model.sigma_T.sigma_f(u), (1,:))
     alpha = quanto_drift(m.factor_alias, m.hjm_model.quanto_model, s, t, X)
     #
     chi = chi_hjm(m)
@@ -208,7 +210,7 @@ function Sigma_T(
     @assert isnothing(X) == !state_dependent_Sigma(m)
     chi = chi_hjm(m)
     # make sure we do not apply correlations twice!
-    f(u) = H_hjm(chi,u,t) .* m.hjm_model.sigma_T.HHfInv .* reshape(m.hjm_model.sigma_T.sigma_f(u), (1,:))
+    f(u) = H_hjm(chi,u,t) .* m.hjm_model.sigma_T.scaling_matrix .* reshape(m.hjm_model.sigma_T.sigma_f(u), (1,:))
     return f
 end
 

src/models/rates/GaussianHjmModel.jl

@@ -1,25 +1,25 @@
 
 """
     struct GaussianHjmModelVolatility
-        HHfInv::AbstractMatrix
+        scaling_matrix::AbstractMatrix
         sigma_f::BackwardFlatVolatility
         DfT::AbstractMatrix
     end
 
 A dedicated matrix-valued volatility term structure for Gaussian HJM Models.
 """
 struct GaussianHjmModelVolatility
-    HHfInv::AbstractMatrix
+    scaling_matrix::AbstractMatrix
     sigma_f::BackwardFlatVolatility
     DfT::AbstractMatrix
 end
 
 function volatility(o::GaussianHjmModelVolatility, u::ModelTime)
-    # return o.HHfInv * (o.DfT .* o.sigma_f(u))  # beware DfT multiplication
+    # return o.scaling_matrix * (o.DfT .* o.sigma_f(u))  # beware DfT multiplication
     σ = o.sigma_f(u)
     d = length(σ)
     return [
-        sum( o.HHfInv[i,k] * o.DfT[k,j] * σ[k] for k = 1:d )
+        sum( o.scaling_matrix[i,k] * o.DfT[k,j] * σ[k] for k = 1:d )
         for i = 1:d, j = 1:d
     ]
 end
@@ -64,6 +64,7 @@ end
         sigma_f::BackwardFlatVolatility,
         correlation_holder::Union{CorrelationHolder, Nothing},
         quanto_model::Union{AssetModel, Nothing},
+        scaling_type::BenchmarkTimesScaling = ForwardRateScaling,
         )
 
 Create a Gausian HJM model.
@@ -75,6 +76,7 @@ function gaussian_hjm_model(
     sigma_f::BackwardFlatVolatility,
     correlation_holder::Union{CorrelationHolder, Nothing},
     quanto_model::Union{AssetModel, Nothing},
+    scaling_type::BenchmarkTimesScaling = ForwardRateScaling,
     )
     # Check inputs
     @assert length(delta()) > 0
@@ -101,8 +103,8 @@ function gaussian_hjm_model(
         DfT = cholesky(Gamma).L
     end
     # prepare vol calculation
-    HHfInv = benchmark_times_scaling(chi(),delta())
-    sigma_T = GaussianHjmModelVolatility(HHfInv, sigma_f, DfT)
+    scaling_matrix = benchmark_times_scaling(chi(), delta(), scaling_type)
+    sigma_T = GaussianHjmModelVolatility(scaling_matrix, sigma_f, DfT)
     # pre-calculate variance 
     d = length(delta())
     y = zeros(d, d, 0)
@@ -216,7 +218,7 @@ function Theta(
     @assert isnothing(X) == !state_dependent_Theta(m)
     y(t) = func_y(m, t)
     # make sure we do not apply correlations twice in quanto adjustment!
-    sigma_T_hyb(u) = m.sigma_T.HHfInv .* reshape(m.sigma_T.sigma_f(u), (1,:))
+    sigma_T_hyb(u) = m.sigma_T.scaling_matrix .* reshape(m.sigma_T.sigma_f(u), (1,:))
     alpha = quanto_drift(m.factor_alias, m.quanto_model, s, t, X)
     return vcat(
         func_Theta_x_integrate_y(m.chi(),y,sigma_T_hyb,alpha,s,t,parameter_grid(m)),
@@ -265,7 +267,7 @@ function Sigma_T(
     )
     @assert isnothing(X) == !state_dependent_Sigma(m)
     # make sure we do not apply correlations twice!
-    sigma_T_hyb(u) = m.sigma_T.HHfInv .* reshape(m.sigma_T.sigma_f(u), (1,:))
+    sigma_T_hyb(u) = m.sigma_T.scaling_matrix .* reshape(m.sigma_T.sigma_f(u), (1,:))
     return func_Sigma_T(m.chi(),sigma_T_hyb,s,t)
 end
 

src/models/rates/SeparableHjmModel.jl

@@ -91,18 +91,84 @@ function G_hjm(m::SeparableHjmModel, s::ModelTime, T::AbstractVector)
     return G_hjm(chi_hjm(m), s, T)
 end
 
+
+"""
+    @enum(
+        BenchmarkTimesScaling,
+        ForwardRateScaling,
+        ZeroRateScaling,
+        DiagonalScaling,
+    )
+
+Specify scaling methods to be used for state variable `x` volatility function.
+
+`ForwardRateScaling` (default) uses forward rates as benchmark rates. For refrence,
+see Andersen/Piterbarg, Interest Rate Modeling, 2010, Prop. 13.3.2.
+
+`ZeroRateScaling` uses continous compounded zero rates as benchmark rates. See
+https://ssrn.com/abstract=4638188 for details.
+
+`DiagonalScaling` uses identity matrix (i.e. no scaling).
 """
-    benchmark_times_scaling(chi::AbstractVector, delta::AbstractVector)
+@enum(
+    BenchmarkTimesScaling,
+    ForwardRateScaling,
+    ZeroRateScaling,
+    DiagonalScaling,
+)
+
+
+"""
+    benchmark_times_scaling_forward_rate(chi::AbstractVector, delta::AbstractVector)
 
 Benchmark times volatility scaling matrix ``H [H^f]^{-1} = [H^f H^{-1}]^{-1}``.
 """
-function benchmark_times_scaling(chi::AbstractVector, delta::AbstractVector)
+function benchmark_times_scaling_forward_rate(chi::AbstractVector, delta::AbstractVector)
     # Hf_H_inv = exp.(-delta * chi')
     Hf_H_inv = [ exp(-chi_ * delta_) for delta_ in delta, chi_ in chi ]  # beware the order of loops!
     HHfInv = inv(Hf_H_inv)
     return HHfInv
 end
 
+
+"""
+    benchmark_times_scaling_zero_rate(chi::AbstractVector, delta::AbstractVector)
+
+Benchmark times volatility scaling matrix based on zero rates.
+"""
+function benchmark_times_scaling_zero_rate(chi::AbstractVector, delta::AbstractVector)
+    A_tau = [ (1.0 - exp(-chi_ * delta_)) / (chi_ * delta_) for delta_ in delta, chi_ in chi ]  # beware the order of loops!
+    return inv(A_tau)
+end
+
+
+"""
+    benchmark_times_scaling(
+        chi::AbstractVector,
+        delta::AbstractVector,
+        scaling_type::BenchmarkTimesScaling = ForwardRate::BenchmarkTimesScaling
+        )
+
+Calculate volatility scaling matrix depending on the scaling type.
+"""
+function benchmark_times_scaling(
+    chi::AbstractVector,
+    delta::AbstractVector,
+    scaling_type::BenchmarkTimesScaling = ForwardRateScaling,
+    )
+    #
+    if scaling_type == ForwardRateScaling
+        return benchmark_times_scaling_forward_rate(chi, delta)
+    elseif scaling_type == ZeroRateScaling
+        return benchmark_times_scaling_zero_rate(chi, delta)
+    elseif scaling_type == DiagonalScaling
+        return Matrix(I, length(delta), length(delta))
+    else
+        error("Unknown BenchmarkTimesScaling type " * string(scaling_type))
+    end
+end
+
+
 """
     func_y(
         y0::AbstractMatrix,

test/unittests/models/futures/markov_future_model.jl

@@ -18,11 +18,23 @@ using Test
         @test DiffFusion.alias(m) == "Std"
         @test DiffFusion.state_alias(m)  == [ "Std_x_1", "Std_x_2", "Std_x_3", ]
         @test DiffFusion.factor_alias(m) == [ "Std_f_1", "Std_f_2", "Std_f_3", ]
+        #
         HHfInv = DiffFusion.benchmark_times_scaling(chi(),delta())
         @test m.hjm_model.sigma_T(1.5) == HHfInv * Diagonal([60., 70., 80.])
         @test m.hjm_model.sigma_T(5.0) == HHfInv * Diagonal([70., 80., 90.])
         @test m.hjm_model.sigma_T(12.0) == HHfInv * Diagonal([80., 90., 90.])
         #
+        m = DiffFusion.markov_future_model("Std",delta,chi,sigma_F,nothing,nothing,DiffFusion.ZeroRateScaling)
+        A_inv = DiffFusion.benchmark_times_scaling(chi(), delta(), DiffFusion.ZeroRateScaling)
+        @test m.hjm_model.sigma_T(1.5) == A_inv * Diagonal([60., 70., 80.])
+        @test m.hjm_model.sigma_T(5.0) == A_inv * Diagonal([70., 80., 90.])
+        @test m.hjm_model.sigma_T(12.0) == A_inv * Diagonal([80., 90., 90.])
+        #
+        m = DiffFusion.markov_future_model("Std",delta,chi,sigma_F,nothing,nothing,DiffFusion.DiagonalScaling)
+        @test m.hjm_model.sigma_T(1.5) == Diagonal([60., 70., 80.])
+        @test m.hjm_model.sigma_T(5.0) == Diagonal([70., 80., 90.])
+        @test m.hjm_model.sigma_T(12.0) == Diagonal([80., 90., 90.])
+        #
         @test_throws AssertionError DiffFusion.markov_future_model("Std", DiffFusion.flat_parameter("", delta()[2:end]), chi, sigma_F, nothing, nothing)
         @test_throws AssertionError DiffFusion.markov_future_model("Std", delta, DiffFusion.flat_parameter("", chi()[2:end]), sigma_F, nothing, nothing)
         @test_throws AssertionError DiffFusion.markov_future_model("Std", DiffFusion.flat_parameter("", reverse(delta())), chi, sigma_F, nothing, nothing)

test/unittests/models/rates/gaussian_hjm_model.jl

@@ -19,10 +19,23 @@ using LinearAlgebra
         @test DiffFusion.alias(m) == "Std"
         @test DiffFusion.state_alias(m) == [ "Std_x_1", "Std_x_2", "Std_x_3", "Std_s" ]
         @test DiffFusion.factor_alias(m) == [ "Std_f_1", "Std_f_2", "Std_f_3" ]
+        #
         HHfInv = DiffFusion.benchmark_times_scaling(chi(),delta())
         @test m.sigma_T(1.5) == HHfInv * Diagonal([60., 70., 80.])
         @test m.sigma_T(5.0) == HHfInv * Diagonal([70., 80., 90.])
         @test m.sigma_T(12.0) == HHfInv * Diagonal([80., 90., 90.])
+        #
+        m = DiffFusion.gaussian_hjm_model("Std",delta,chi,sigma_f,nothing,nothing,DiffFusion.ZeroRateScaling)
+        A_inv = DiffFusion.benchmark_times_scaling(chi(), delta(), DiffFusion.ZeroRateScaling)
+        @test m.sigma_T(1.5) == A_inv * Diagonal([60., 70., 80.])
+        @test m.sigma_T(5.0) == A_inv * Diagonal([70., 80., 90.])
+        @test m.sigma_T(12.0) == A_inv * Diagonal([80., 90., 90.])
+        #
+        m = DiffFusion.gaussian_hjm_model("Std",delta,chi,sigma_f,nothing,nothing,DiffFusion.DiagonalScaling)
+        @test m.sigma_T(1.5) == Diagonal([60., 70., 80.])
+        @test m.sigma_T(5.0) == Diagonal([70., 80., 90.])
+        @test m.sigma_T(12.0) == Diagonal([80., 90., 90.])
+        #
         @test DiffFusion.correlation_holder(m) == nothing
         #
         @test_throws AssertionError DiffFusion.gaussian_hjm_model("Std", DiffFusion.flat_parameter("", delta()[2:end]), chi, sigma_f, nothing, nothing)

test/unittests/models/rates/separable_hjm_model.jl

@@ -39,7 +39,10 @@ using LinearAlgebra
             chi[1]*delta[2] chi[2]*delta[2] chi[3]*delta[2];
             chi[1]*delta[3] chi[2]*delta[3] chi[3]*delta[3];
         ]
-        @test DiffFusion.benchmark_times_scaling(chi,delta) == inv(exp.(-chi_delta))
+        @test DiffFusion.benchmark_times_scaling(chi, delta) == inv(exp.(-chi_delta))
+        @test DiffFusion.benchmark_times_scaling(chi, delta, DiffFusion.ForwardRateScaling) == inv(exp.(-chi_delta))
+        @test DiffFusion.benchmark_times_scaling(chi, delta, DiffFusion.ZeroRateScaling) == inv((1.0 .- exp.(-chi_delta)) ./ chi_delta)
+        @test DiffFusion.benchmark_times_scaling(chi, delta, DiffFusion.DiagonalScaling) == Matrix(I, 3, 3)
     end
 
     @testset "Auxilliary state variable/variance calculation." begin