no need for storage

matthieugomez · Aug 30, 2021 · f74e9ac · f74e9ac · matthieugomez · Aug 30, 2021
1 parent f625f8e
commit f74e9ac
Show file tree

Hide file tree

Showing 2 changed files with 24 additions and 35 deletions.
diff --git a/src/AdditiveFunctional.jl b/src/AdditiveFunctional.jl
@@ -23,11 +23,9 @@ end
 compute the tail index of the stationary distribution of e^{m}, i.e.
 ζ such that cgf(m)(ζ) = δ
 """
-function tail_index(m::AdditiveFunctional; δ = 0.0, verbose = false, r0 = ones(length(m.X.x)), xatol = 1e-4, kwargs...)
-    r0 = deepcopy(r0)
+function tail_index(m::AdditiveFunctional; δ = 0, verbose = false, r0 = Ones(length(m.X.x)), xatol = 1e-4, kwargs...)
     fzero((1e-5, 1e3); xatol = xatol, kwargs...) do ξ
-        η, f = cgf(m; r0 = r0)(ξ)
-        copyto!(r0, f)
+        η, r0 = cgf(m; r0 = r0)(ξ)
         verbose && @show (:LR, ξ, η)
         return η - δ
     end
@@ -58,13 +56,13 @@ mutable struct AdditiveFunctionalDiffusion <: AdditiveFunctional
     μm::AbstractVector{<:Number}
     σm::AbstractVector{<:Number}
     ρ::Number
-    𝕋::Tridiagonal
 end
 
 function AdditiveFunctionalDiffusion(X::DiffusionProcess, μm::AbstractVector{<:Number}, σm::AbstractVector{<:Number}; ρ::Number = 0.0)
-    AdditiveFunctionalDiffusion(X, μm, σm, ρ, deepcopy(X.𝕋))
+    length(X.x) == length(μm) == length(σm) || throw(ArgumentError("Vector for grid, drift, and volatility should have the same size"))
+    AdditiveFunctionalDiffusion(X, μm, σm, ρ)
 end
 
 function generator(M::AdditiveFunctionalDiffusion)
-    ξ -> Diagonal(ξ .* M.μm .+ 0.5 * ξ^2 .* M.σm.^2) + generator!(M.𝕋, M.X.x, M.X.μx .+ ξ .* M.ρ .* M.σm .* M.X.σx, M.X.σx)
+    ξ -> generator(M.X.x, ξ .* M.μm .+ 0.5 * ξ^2 .* M.σm.^2, M.X.μx .+ ξ .* M.ρ .* M.σm .* M.X.σx, M.X.σx)
 end
diff --git a/src/MarkovProcess.jl b/src/MarkovProcess.jl
@@ -29,26 +29,24 @@ mutable struct DiffusionProcess <: MarkovProcess
     x::AbstractVector{<:Real}
     μx::AbstractVector{<:Real}
     σx::AbstractVector{<:Real}
-    𝕋::Tridiagonal
+    function DiffusionProcess(x::AbstractVector{<:Real}, μx::AbstractVector{<:Real}, σx::AbstractVector{<:Real})
+        length(x) == length(μx) == length(σx) || throw(ArgumentError("Vector for grid, drift, and volatility should have the same size"))
+        new(x, μx, σx)
+    end
 end
 
-function DiffusionProcess(x::AbstractVector{<:Real}, μx::AbstractVector{<:Real}, σx::AbstractVector{<:Real})
-    length(x) == length(μx) || error("Vector for grid, drift, and volatility should have the same size")
-    length(μx) == length(σx) || error("Vector for grid, drift, and volatility should have the same size")
-    n = length(x)
-    𝕋 = Tridiagonal(zeros(n-1), zeros(n), zeros(n-1))
-    generator!(𝕋, x, μx, σx)
-    DiffusionProcess(x, μx, σx, 𝕋)
-end
 
-generator(X::DiffusionProcess) = X.𝕋
 state_space(X::DiffusionProcess) = X.x
 
-# Compute the generator 
-function generator!(𝕋, x, μx::AbstractVector, σx::AbstractVector)
+function generator(X::DiffusionProcess) 
+    generator(X.x, Zeros(length(X.x)), X.μx, X.σx)
+end
+
+# create operator associated with f ⭌ v * f + μx * ∂f + 0.5 * σx^2 * ∂^2f
+function generator(x::AbstractVector, v::AbstractVector, μx::AbstractVector, σx::AbstractVector)
     # The key is that sum of each row = 0.0 and off diagonals are positive
     n = length(x)
-    fill!(𝕋, 0)
+    𝕋 = Tridiagonal(zeros(n-1), zeros(n), zeros(n-1))
     @inbounds for i in 1:n
         Δxp =x[min(i, n-1)+1] - x[min(i, n-1)]
         Δxm = x[max(i-1, 1) + 1] - x[max(i-1, 1)]
@@ -68,15 +66,17 @@ function generator!(𝕋, x, μx::AbstractVector, σx::AbstractVector)
     # ensure machine precision
     c = sum(𝕋, dims = 2)
     for i in 1:n
-        𝕋[i, i] -= c[i]
+        𝕋[i, i] += v[i] - c[i]
     end
     return 𝕋
 end
 
+# create operator associated with f ⭌ ∂f using upwinding w.r.t. μx
+function ∂(X::DiffusionProcess)
+    Diagonal(X.μx) \ generator(X.x, Zeros(length(X.x)), X.μx, Zeros(length(X.x)))
+end
 
-
-
-# Special cases.
+# Special Diffusion Processes
 # it's important to take low p to have the right tail index of Additive functional (see tests)
 function OrnsteinUhlenbeck(; xbar = 0.0, κ = 0.1, σ = 1.0, p = 1e-10, length = 100, 
     xmin = quantile(Normal(xbar, σ / sqrt(2 * κ)), p), xmax = quantile(Normal(xbar, σ / sqrt(2 * κ)), 1 - p))
@@ -85,21 +85,12 @@ function OrnsteinUhlenbeck(; xbar = 0.0, κ = 0.1, σ = 1.0, p = 1e-10, length =
     else
         x = range(xmin, stop = xmax, length = length)
     end
-    μx = κ .* (xbar .- x)
-    σx = σ .* Ones(Base.length(x))
-    DiffusionProcess(x, μx, σx)
+    DiffusionProcess(x, κ .* (xbar .- x), σ * Ones(Base.length(x)))
 end
 
 function CoxIngersollRoss(; xbar = 0.1, κ = 0.1, σ = 1.0, p = 1e-10, length = 100, α = 2 * κ * xbar / σ^2, β = σ^2 / (2 * κ), xmin = quantile(Gamma(α, β), p), xmax = quantile(Gamma(α, β), 1 - p), pow = 2)
     # check 0 is not attainable
     @assert (2 * κ * xbar) / σ^2 > 1
     x = range(xmin^(1/pow), stop = xmax^(1/pow), length = length).^pow
-    μx = κ .* (xbar .- x)
-    σx = σ .* sqrt.(x)
-    DiffusionProcess(x, μx, σx)
-end
-
-
-function ∂(X::DiffusionProcess)
-    Diagonal(X.μx) \ generator!(deepcopy(X.𝕋), X.x, X.μx, Zeros(length(X.x)))
+    DiffusionProcess(x, κ .* (xbar .- x), σ .* sqrt.(x))
 end