T is back to being generator

matthieugomez · Mar 17, 2020 · d4e9969 · d4e9969
1 parent cce7236
commit d4e9969
Showing 1 changed file with 22 additions and 13 deletions.
diff --git a/src/diffusions.jl b/src/diffusions.jl
@@ -1,8 +1,10 @@
 #========================================================================================
 
-Compute the operator
+Compute the operator 
 Tf = v_0 * f + v1 * ∂(f) + v2 * ∂∂(f)
-Note that
+
+
+Mote that
 T'g = v_0 * g - ∂(v1 * g) + ∂∂(v2 * g)
 
 ========================================================================================#
@@ -22,24 +24,31 @@ function make_Δ(x)
     return x, 1 ./ Δx, 1 ./ Δxm, 1 ./ Δxp
 end
 
-function build_diffusion!(T, Δ, v0::AbstractVector, v1::AbstractVector, v2::AbstractVector)
+function operator_second_order(x::AbstractVector, v0::AbstractVector, v1::AbstractVector, v2::AbstractVector)
+    n = length(x)
+    T = Tridiagonal(zeros(n-1), zeros(n), zeros(n-1))
+    Δ = make_Δ(x)
+    operator_second_order!(T, Δ, v0, v1, v2)
+end
+
+function operator_second_order!(T, Δ, v0::AbstractVector, v1::AbstractVector, v2::AbstractVector)
     # The key is that sum of each row = 0.0 and off diagonals are positive (singular M-matrix)
     x, invΔx, invΔxm, invΔxp = Δ
     n = length(x)
     fill!(T, 0)
     @inbounds for i in 1:n
         if v1[i] >= 0
-            T[min(i + 1, n), i] += v1[i] * invΔxp[i]
+            T[i, min(i + 1, n)] += v1[i] * invΔxp[i]
             T[i, i] -= v1[i] * invΔxp[i]
         else
             T[i, i] += v1[i] * invΔxm[i]
-            T[max(i - 1, 1), i] -= v1[i] * invΔxm[i]
+            T[i, max(i - 1, 1)] -= v1[i] * invΔxm[i]
         end
-        T[max(i - 1, 1), i] += v2[i] * invΔxm[i] * invΔx[i]
+        T[i, max(i - 1, 1)] += v2[i] * invΔxm[i] * invΔx[i]
         T[i, i] -= v2[i] * 2 * invΔxm[i] * invΔxp[i]
-        T[min(i + 1, n), i] += v2[i] * invΔxp[i] * invΔx[i]
+        T[i, min(i + 1, n)] += v2[i] * invΔxp[i] * invΔx[i]
     end
-    c = sum(T, dims = 1)
+    c = sum(adjoint(T), dims = 1)
     for i in 1:n
         T[i, i] += v0[i] - c[i]
     end
@@ -67,7 +76,7 @@ function DiffusionProcess(x::AbstractVector{<:Real}, μx::AbstractVector{<:Real}
     n = length(x)
     T = Tridiagonal(zeros(n-1), zeros(n), zeros(n-1))
     Δ = make_Δ(x)
-    build_diffusion!(T, Δ, Zeros(length(x)), μx, 0.5 * σx.^2)
+    operator_second_order!(T, Δ, Zeros(length(x)), μx, 0.5 * σx.^2)
     DiffusionProcess(x, μx, σx, T, Δ)
 end
 
@@ -91,12 +100,12 @@ function CoxIngersollRoss(; xbar = 0.1, κ = 0.1, σ = 1.0, p = 1e-10, length =
     DiffusionProcess(x, μx, σx)
 end
 
-generator(X::DiffusionProcess) = adjoint(X.T)
+generator(X::DiffusionProcess) = X.T
 state_space(X::DiffusionProcess) = X.x
 
 function ∂(X::DiffusionProcess)
-    T = build_diffusion!(deepcopy(X.T), X.Δ, Zeros(length(X.x)), Ones(length(X.x)), Zeros(length(X.x)))
-    adjoint(T)
+    T = deepcopy(X.T)
+    operator_second_order!(T, X.Δ, Zeros(length(X.x)), Ones(length(X.x)), Zeros(length(X.x)))
 end
 
 #========================================================================================
@@ -123,5 +132,5 @@ function AdditiveFunctionalDiffusion(X::DiffusionProcess, μm::AbstractVector{<:
 end
 
 function generator(M::AdditiveFunctionalDiffusion)
-    ξ -> build_diffusion!(M.T, M.X.Δ, ξ .* M.μm .+ 0.5 * ξ^2 .* M.σm.^2 .- M.δ,  M.X.μx .+ ξ .* M.ρ .* M.σm .* M.X.σx, 0.5 * M.X.σx.^2)'
+    ξ -> operator_second_order!(M.T, M.X.Δ, ξ .* M.μm .+ 0.5 * ξ^2 .* M.σm.^2 .- M.δ,  M.X.μx .+ ξ .* M.ρ .* M.σm .* M.X.σx, 0.5 * M.X.σx.^2)
 end