Correlated Generators (#39)

Project.toml

@@ -4,9 +4,12 @@ authors = ["Alec Loudenback <alecloudenback@gmail.com> and contributors"]
 version = "0.3.8"
 
 [deps]
+Copulas = "ae264745-0b69-425e-9d9d-cf662c5eec93"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 Yields = "d7e99b2f-e7f3-4d9e-9f01-2338fc023ad3"
 
 [compat]

README.md

@@ -187,6 +187,35 @@ p
 
 ![BSM Paths](https://user-images.githubusercontent.com/711879/180128072-fe08d285-0edc-4707-a89e-8d14fef23d2a.png)
 
+## Correlated Scenarios
+
+Combined with `using Copulas`, you can create correlated scenarios with a given copula. See `?Correlated` for the docstring on creating a correlated set of scenario generators.
+
+### Example
+
+Create two equity paths that are 90% correlated:
+
+```
+using EconomicScenarioGenerators, Copulas
+
+m = BlackScholesMerton(0.01,0.02,.15,100.)
+s = ScenarioGenerator(
+                      1,  # timestep
+                      30, # projection horizon
+                      m,  # model
+                  )
+
+ss = [s,s] # these don't have to be the exact same, but do need same shape
+g = GaussianCopula([1. 0.9; 0.9 1.]) # 90% correlated
+c = Correlated(ss,g)
+
+collect(c)
+
+using Plots
+plot(collect(c))
+```
+![plot_2](https://user-images.githubusercontent.com/711879/183233379-c9dda6ba-c945-4e69-937d-aa7835198f5e.svg)
+
 ## Benchmarks
 
 Generating 10,000 scenarios of daily timesteps for 1 year (252 business days) with a Black-Scholes-Merton model:

docs/src/index.md

@@ -193,6 +193,35 @@ p
 
 ![BSM Paths](https://user-images.githubusercontent.com/711879/180128072-fe08d285-0edc-4707-a89e-8d14fef23d2a.png)
 
+## Correlated Scenarios
+
+Combined with `using Copulas`, you can create correlated scenarios with a given copula. See `?Correlated` for the docstring on creating a correlated set of scenario generators.
+
+### Example
+
+Create two equity paths that are 90% correlated:
+
+```
+using EconomicScenarioGenerators, Copulas
+
+m = BlackScholesMerton(0.01,0.02,.15,100.)
+s = ScenarioGenerator(
+                      1,  # timestep
+                      30, # projection horizon
+                      m,  # model
+                  )
+
+ss = [s,s] # these don't have to be the exact same, but do need same shape
+g = GaussianCopula([1. 0.9; 0.9 1.]) # 90% correlated
+c = Correlated(ss,g)
+
+collect(c)
+
+using Plots
+plot(collect(c))
+```
+![plot_2](https://user-images.githubusercontent.com/711879/183233379-c9dda6ba-c945-4e69-937d-aa7835198f5e.svg)
+
 ## Benchmarks
 
 Generating 10,000 scenarios of daily timesteps for 1 year (252 business days) with a Black-Scholes-Merton model:

src/EconomicScenarioGenerators.jl

@@ -4,6 +4,9 @@ import ForwardDiff
 import Yields
 import IterTools
 using Random
+using Copulas
+using Distributions
+using SplitApplyCombine
 
 abstract type EconomicModel end
 
@@ -20,7 +23,29 @@ end
 include("interest.jl")
 include("equity.jl")
 
-struct ScenarioGenerator{N<:Real,T,R<:AbstractRNG}
+abstract type AbstractScenarioGenerator end
+
+
+"""
+    ScenarioGenerator(timestep::N,endtime::N,model,RNG::AbstractRNG) where {N<:Real}
+
+A `ScenarioGenerator` is an iterator which yields the time series of the model. It takes the parameters of the scenario generation such as `timestep` and `endtime` (the time horizon). `model` is any EconomicModel from the `EconomicScenarioGenerators` package.
+
+A `ScenarioGenerator` can be `iterate`d or `collect`ed. 
+
+# Examples
+
+```julia
+m = Vasicek(0.136,0.0168,0.0119,Continuous(0.01)) # a, b, σ, initial Rate
+s = ScenarioGenerator(
+        1,  # timestep
+        30, # projection horizon
+        m,  # model
+    )
+```
+
+"""
+struct ScenarioGenerator{N<:Real,T,R<:AbstractRNG} <: AbstractScenarioGenerator
     timestep::N
     endtime::N
     model::T
@@ -45,10 +70,11 @@ function Base.iterate(sg::ScenarioGenerator{N,T,R},state) where {N,T<:EconomicMo
     if (state.time > sg.endtime) || (state.time ≈ sg.endtime)
         return nothing
     else
-        new_rate = nextrate(sg.RNG,sg.model,state.value,state.time,sg.timestep)
+        variate = rand(sg.RNG) # a quantile
+        new_value = nextvalue(sg.model,state.value,state.time,sg.timestep,variate)
         state = (
             time = state.time + sg.timestep,
-            value = new_rate
+            value = new_value
         )
         return (state.value, state)
     end
@@ -69,9 +95,95 @@ end
 
 Base.eltype(::Type{ScenarioGenerator{N,T,R}}) where {N,T,R} = __outputtype(T)
 
+
+"""
+    Correlated(v::Vector{ScenarioGenerator},copula,RNG::AbstractRNG)
+
+An iterator which uses the given copula to correlate the outcomes of the underling ScenarioGenerators. The copula returns the sampled CDF which the individual models will interpret according to their own logic (e.g. the random variate for the BlackScholesMerton model assumes a normal variate within the diffusion process).
+
+The `time` and `timestep` of the underlying ScenarioGenerators are asserted to be the same upon construction.
+
+A `Correlated` can be `iterate`d or `collect`ed. It will iterate through each sub-scenario generator and yield the time series of each (as opposed to iterating through each timestep for all models at each step).
+
+# Examples
+
+```julia
+using EconomicScenarioGenerators, Copulas
+
+m = BlackScholesMerton(0.01,0.02,.15,100.)
+s = ScenarioGenerator(
+                      1,  # timestep
+                      30, # projection horizon
+                      m,  # model
+                  )
+
+ss = [s,s] # these don't have to be the exact same
+g = GaussianCopula([1. 0.9; 0.9 1.])
+c = Correlated(ss,g)
+
+# collect a vector of two correlated paths.
+collect(c)
+```
+
+"""
+struct Correlated{T,U,R} <: AbstractScenarioGenerator
+    sg::Vector{T}
+    copula::U
+    RNG::R
+    function Correlated(generators::Vector{T},copula::U,RNG::R=Random.GLOBAL_RNG) where {T<:ScenarioGenerator,U,R<:AbstractRNG}
+        fst = first(generators)
+        @assert all(fst.timestep == g.timestep for g in generators) "All component generators must have the same `timestep`."
+        @assert all(fst.endtime == g.endtime for g in generators) "All component generators must have the same `endpoint`."
+        new{T,U,R}(generators,copula,RNG)
+    end
+end
+Base.Broadcast.broadcastable(x::T) where {T<:Correlated} = Ref(x)
+
+function Base.iterate(sgc::Correlated) 
+    n = 1
+    sg = sgc.sg[n]
+    variates = rand(sgc.RNG,sgc.copula,length(sg)) # CDF
+    times = sg.timestep:sg.timestep:(sg.endtime+sg.timestep)
+    scenariovalue = initial_value(sg.model,first(times))
+    values = map(enumerate(times)) do (i,t)
+        scenariovalue = nextvalue(sg.model,scenariovalue,t,sg.timestep,variates[n,i])
+        scenariovalue 
+    end
+
+    state = (
+        variates = variates,
+        n = 2,
+        )
+
+    return (values,state)
+end
+
+function Base.iterate(sgc::Correlated,state)
+    if state.n > length(sgc.sg)
+        return nothing
+    else
+        n = state.n
+        sg = sgc.sg[n]
+        times = sg.timestep:sg.timestep:(sg.endtime+sg.timestep)
+        scenariovalue = initial_value(sg.model,first(times))
+        values = map(enumerate(times)) do (i,t)
+            scenariovalue = nextvalue(sg.model,scenariovalue,t,sg.timestep,state.variates[n,i])
+            scenariovalue 
+        end
+        state = (
+            variates = state.variates,
+            n = n+1,
+            )
+        return (values, state)
+    end
+end
+
+Base.length(sgc::Correlated) = length(sgc.sg)
+Base.eltype(::Type{Correlated{N,T,R}}) where {N,T,R} = Vector{eltype(N)}
+
 include("Yields.jl")
 export Vasicek, CoxIngersollRoss, HullWhite,
         BlackScholesMerton,ConstantElasticityofVariance,
-        ScenarioGenerator, YieldCurve
+        ScenarioGenerator, YieldCurve, Correlated
 
 end

src/Yields.jl

@@ -4,6 +4,16 @@ function YieldCurve(sg::ScenarioGenerator{N,T,R}) where {N,T<:InterestRateModel,
 
 end
 
+function YieldCurve(C::T) where {T<:Correlated}
+    timestep = first(C.sg).timestep
+    _, times = __zeros_times(first(C.sg))
+    _zeros(v) = cumsum(v .* timestep) ./ times
+    rates = collect(C)
+    zs = _zeros.(rates)
+    map(r->Yields.Zero(r,times),zs)
+
+end
+
 function __disc_rate_to_fwd(rate,time)
     log(rate) / -time
 end

src/equity.jl

@@ -35,13 +35,17 @@ struct BlackScholesMerton{T,U,V} <:EquityModel
     σ::V # roughly equivalent to the volatility in the usual lognormal model multiplied by F^{1-β}_{0}
     initial::Float64
 end
-function nextrate(RNG,M::BlackScholesMerton,prior,time,timestep)
+
+
+
+function nextvalue(M::BlackScholesMerton,prior,time,timestep,variate)
+    variate = quantile(Normal(),variate)
     r, q, σ = M.r, M.q, M.σ
     # Hull Options, Futures, & Other Derivatives, 10th ed., pg 470
     return (
         prior *
         exp((r- q - σ^2 / 2) *  timestep + 
-         σ * sqrt(timestep) * randn(RNG))
+         σ * sqrt(timestep) * variate)
     )
 end
 

src/interest.jl

@@ -67,8 +67,9 @@ struct Vasicek{T<:Yields.Rate} <: ShortRateModel
     initial::T # 0.01
 end
 
-function nextrate(RNG,M::Vasicek{T},prior,time,timestep) where T
-    prior + M.a * (M.b - prior) * timestep + M.σ * sqrt(timestep) * randn(RNG)
+function nextvalue(M::Vasicek{T},prior,time,timestep,variate) where T
+    variate = quantile(Normal(),variate)
+    prior + M.a * (M.b - prior) * timestep + M.σ * sqrt(timestep) * variate
 end
 
 """
@@ -139,8 +140,9 @@ struct CoxIngersollRoss{T<:Yields.Rate} <: ShortRateModel
     initial::T # 0.01
 end
 
-function nextrate(RNG,M::CoxIngersollRoss{T},prior,time,timestep) where {T}
-    prior + M.a * (M.b - prior) * timestep + M.σ * sqrt(Yields.rate(prior)) * randn(RNG)
+function nextvalue(M::CoxIngersollRoss{T},prior,time,timestep,variate) where {T}
+    variate = quantile(Normal(),variate)
+    prior + M.a * (M.b - prior) * timestep + M.σ * sqrt(Yields.rate(prior)) * variate
 end
 
 """
@@ -158,10 +160,11 @@ end
 # See Yields.jl for HullWhite with a YieldCurve defining theta
 
 # how would HullWhite be constructed if not giving it a curve?
-function nextrate(RNG,M::HullWhite{T},prior,time,timestep) where {T}
+function nextvalue(M::HullWhite{T},prior,time,timestep,variate) where {T}
+    variate = quantile(Normal(),variate)
     θ_t = θ(M,time+timestep,timestep)
     # https://quantpie.co.uk/srm/hull_white_sr.php
-    prior + (θ_t - M.a * prior) * timestep + M.σ * √(timestep) * randn(RNG)
+    prior + (θ_t - M.a * prior) * timestep + M.σ * √(timestep) * variate
 end
 
 function θ(M::HullWhite{T},time,timestep) where {T<:Yields.AbstractYield}

test/Correlated.jl

@@ -0,0 +1,83 @@
+#TODO test that non-congruent component generators are not allowed
+@testset "Correlated" begin
+    using EconomicScenarioGenerators, Copulas
+
+    @testset "Equity" begin
+        @testset "BSM" begin
+            m = BlackScholesMerton(0.01, 0.02, 0.15, 100.0)
+            s = ScenarioGenerator(
+                1,  # timestep
+                30, # projection horizon
+                m,  # model
+            )
+
+            ss = [s, s] # these don't have to be the exact same
+            g = GaussianCopula([1.0 0.9; 0.9 1.0])
+            c = Correlated(ss, g, StableRNG(1))
+
+            x = collect(c)
+
+            ρ = cor(hcat(ratio.(x)...))
+            @test ρ[1, 2] ≈ 0.9 atol = 0.03
+        end
+    end
+
+    @testset "Interest" begin
+        @testset "HullWhite" begin
+            rates = [0.01, 0.01, 0.03, 0.05, 0.07, 0.16, 0.35, 0.92, 1.40, 1.74, 2.31, 2.41] ./ 100
+            mats = [1 / 12, 2 / 12, 3 / 12, 6 / 12, 1, 2, 3, 5, 7, 10, 20, 30]
+            c = Yields.CMT(rates, mats)
+
+            m = HullWhite(0.1, 0.01, c)
+
+            s = ScenarioGenerator(
+                0.1,                              # timestep
+                30.0,                             # projection horizon
+                m,
+                StableRNG(1)
+            )
+
+            g = GaussianCopula([1.0 0.9; 0.9 1.0])
+            c = Correlated([s, s], g, StableRNG(1))
+
+            x = collect(c)
+
+            ρ = cor(hcat(diff.(map(y->Yields.rate.(y),x))...))
+            @test ρ[1, 2] ≈ 0.9 atol = 0.03
+        end
+
+        @testset "Vasicek" begin
+            m = Vasicek(0.136, 0.0168, 0.0119, Yields.Continuous(0.01))
+            s = ScenarioGenerator(
+                .1,                              # timestep
+                30.0,                             # projection horizon
+                m
+            )
+            g = GaussianCopula([1.0 0.9; 0.9 1.0])
+            c = Correlated([s, s], g, StableRNG(1))
+
+            x = collect(c)
+
+            ρ = cor(hcat(diff.(map(y->Yields.rate.(y),x))...))
+            @test ρ[1, 2] ≈ 0.9 atol = 0.03
+        end
+        @testset "CIR" begin
+            m = CoxIngersollRoss(0.136, 0.0168, 0.0119, Yields.Continuous(0.01))
+            s = ScenarioGenerator(
+                .1,                              # timestep
+                30.0,                             # projection horizon
+                m
+            )
+            g = GaussianCopula([1.0 0.9; 0.9 1.0])
+            c = Correlated([s, s], g, StableRNG(1))
+
+            x = collect(c)
+
+            ρ = cor(hcat(diff.(map(y->Yields.rate.(y),x))...))
+            @test ρ[1, 2] ≈ 0.9 atol = 0.03
+
+
+            @test first(YieldCurve.(c)) isa Yields.AbstractYield
+        end
+    end
+end