Simulating CIR model with Differential Equation Package

General Usage
1

I am trying to solve CIR model DiffEq.but get the following error.

Please help me resolving this issue.

# dx=a(b−x)dt+σdWt Vaiscek Model
# dx=a(b−x)dt+σ*sqrt(x)*dWt CIR model

α=0.2;σ=0.2;u₀=0.04;β=0.041
f(u,p,t) = α*(β-u)
g(u,p,t) = σ*sqrt(u)
dt = 1//250 #1//2^(5) # in fractional form
tspan = (0.0, 10.0)
prob = SDEProblem(f,g,u₀,tspan)
sol = solve(prob,EM(),dt=dt)
plot(sol,legend=false)

DomainError with -7.415721101312867e-6:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).

Stacktrace:
 [1] throw_complex_domainerror(::Symbol, ::Float64) at .\math.jl:31
 [2] sqrt at .\math.jl:492 [inlined]
 [3] g(::Float64, ::DiffEqBase.NullParameters, ::Float64) at .\In[3]:9
 [4] perform_step!(::StochasticDiffEq.SDEIntegrator{EM{true},false,Float64,Float64,Float64,DiffEqBase.NullParameters,Float64,Float64,Float64,NoiseProcess{Float64,1,Float64,Float64,Nothing,Nothing,typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE),false,DataStructures.Stack{Tuple{Float64,Float64,Nothing}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Nothing},false},RSWM{:RSwM1,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},Float64,RODESolution{Float64,1,Array{Float64,1},Nothing,Nothing,Array{Float64,1},NoiseProcess{Float64,1,Float64,Float64,Nothing,Nothing,typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE),false,DataStructures.Stack{Tuple{Float64,Float64,Nothing}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Nothing},false},RSWM{:RSwM1,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},SDEProblem{Float64,Tuple{Float64,Float64},false,DiffEqBase.NullParameters,Nothing,SDEFunction{false,typeof(f),typeof(g),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g),Nothing,Nothing},EM{true},StochasticDiffEq.LinearInterpolationData{Array{Float64,1},Array{Float64,1}},DiffEqBase.DEStats},StochasticDiffEq.EMConstantCache,SDEFunction{false,typeof(f),typeof(g),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g),StochasticDiffEq.SDEOptions{Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Float64,Float64,Float64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Nothing,Float64,Nothing}, ::StochasticDiffEq.EMConstantCache, ::Function) at C:\Users\Himanshi\.julia\packages\StochasticDiffEq\e0qFK\src\perform_step\low_order.jl:13
 [5] perform_step! at C:\Users\Himanshi\.julia\packages\StochasticDiffEq\e0qFK\src\perform_step\low_order.jl:2 [inlined]
 [6] solve!(::StochasticDiffEq.SDEIntegrator{EM{true},false,Float64,Float64,Float64,DiffEqBase.NullParameters,Float64,Float64,Float64,NoiseProcess{Float64,1,Float64,Float64,Nothing,Nothing,typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE),false,DataStructures.Stack{Tuple{Float64,Float64,Nothing}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Nothing},false},RSWM{:RSwM1,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},Float64,RODESolution{Float64,1,Array{Float64,1},Nothing,Nothing,Array{Float64,1},NoiseProcess{Float64,1,Float64,Float64,Nothing,Nothing,typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE),false,DataStructures.Stack{Tuple{Float64,Float64,Nothing}},ResettableStacks.ResettableStack{Tuple{Float64,Float64,Nothing},false},RSWM{:RSwM1,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},SDEProblem{Float64,Tuple{Float64,Float64},false,DiffEqBase.NullParameters,Nothing,SDEFunction{false,typeof(f),typeof(g),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g),Nothing,Nothing},EM{true},StochasticDiffEq.LinearInterpolationData{Array{Float64,1},Array{Float64,1}},DiffEqBase.DEStats},StochasticDiffEq.EMConstantCache,SDEFunction{false,typeof(f),typeof(g),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g),StochasticDiffEq.SDEOptions{Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Float64,Float64,Float64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Nothing,Float64,Nothing}) at C:\Users\Himanshi\.julia\packages\StochasticDiffEq\e0qFK\src\solve.jl:401
 [7] #__solve#40 at C:\Users\Himanshi\.julia\packages\StochasticDiffEq\e0qFK\src\solve.jl:7 [inlined]
 [8] #__solve at .\none:0 [inlined] (repeats 5 times)
 [9] #solve#386(::Base.Iterators.Pairs{Symbol,Rational{Int64},Tuple{Symbol},NamedTuple{(:dt,),Tuple{Rational{Int64}}}}, ::Function, ::SDEProblem{Float64,Tuple{Float64,Float64},false,DiffEqBase.NullParameters,Nothing,SDEFunction{false,typeof(f),typeof(g),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g),Nothing,Nothing}, ::EM{true}) at C:\Users\Himanshi\.julia\packages\DiffEqBase\pqp0B\src\solve.jl:39
 [10] (::getfield(DiffEqBase, Symbol("#kw##solve")))(::NamedTuple{(:dt,),Tuple{Rational{Int64}}}, ::typeof(solve), ::SDEProblem{Float64,Tuple{Float64,Float64},false,DiffEqBase.NullParameters,Nothing,SDEFunction{false,typeof(f),typeof(g),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g),Nothing,Nothing}, ::EM{true}) at .\none:0
 [11] top-level scope at In[3]:13
2

EM does not preserve positivity. As you can see, the CIR cannot cross zero because the drift forbids it. You can either use a smaller dt or use Ito formula to solve log(x)

3

This problem can use most methods. Why choose EM here? Likely a higher order adaptive method would better preserve positivity if the real equation is positive

4

I experimented with many methods listed in stochastic solver page in Diffeq documentation but none seem to work particularly when i try to create an ensemble. Even with a smaller dt, it does not work always. I have not tried log(x) version as I need to figure out how to accomplish this.

5

I don’t think an assumption that it has to be positive is appropriate for that equation. Are you sure that the equation you’re solving actually has a positive solution with probability 1? With a mean around 0.041 a standard deviation much larger than 1, first estimates say it has a very high probability that the true solution is negative often negative. Additionally, sqrt(x) noise decreases too slow for it to not hit zero. Remember that u' = -sqrt(u) is an ODE that hits zero in finite time!

Given that it will surely hit zero, sqrt(abs(u)) seems more appropriate then. Or your model parameter might not be what you think they are.

6

I didn’t get you when you say stddev is much larger than 1. The stddev I am using is 0.2 which translates to variance of 0.04. I agree that due to large stddev it may quickly decrease to zero or turn negative. So I tried with stddev of 0.05. It still does not work unless I use sqrt(abs(u)).

I replicated the problem in excel using sqrt(u) and it does not throw any negative there. The equation in excel is approximated as

r(t)=r(t-1)+a(b-r(t-1))*dt+sigma*sqrt(r(t-1))*sqrt(dt)*normsinv(rand())

I am not an expert in using DiffEq tool so may be missing something here.

Update
If I change sqrt(u) to sqrt(u*dt) in the equation even with sigma=0.2, it seems to work.

7

That’s not what I said. I said the standard deviation is larger than the mean, much larger, and so it’s going to hit zero with essentially probability 1.

That’s not your original equation.

You can very easily just run the loop and see that the SDE will go negative.

α=0.2;σ=0.2;u₀=0.04;β=0.041
f(u,p,t) = α*(β-u)
g(u,p,t) = σ*sqrt(u)
dt = 1//250 #1//2^(5) # in fractional form
tspan = (0.0, 10.0)
u = 0.04
for i in 1:10000
  global u
  u = u + dt*f(u,nothing,nothing) + g(u,nothing,nothing)*sqrt(dt)*randn()
  if u < 0
    @show i,u
  end
end

(i, u) = (1520, -5.101939115412398e-5)
(i, u) = (392, -0.00012766187472027155)
(i, u) = (1323, -0.00022494655749871805)
(i, u) = (979, -0.00011945721626705025)
8

This seems to be a pretty well-known result:

http://janroman.dhis.org/finance/IR/CIR%20advanced.pdf

Positivity is guaranteed when

2α*β >= σ^2 #false

The second link there mentions that when the condition is false, you should expect that the solution hits zero infinitely many times as t -> infty with probability 1. So given that the solution should, if calculated correctly, go negative infinitely many times, I would double check that Excel sheet and make sure you translated the same model.

1 Like
9

Thanks. The condition of positivity as mentioned in the same link is 2ab>=sigma^2.

I tested the solution with sqrt(u*dt) as mentioned in the link under Calibration header as I presumed (my ignorance of noise term in DiffEq package) that sqrt(dt) may have to be added to the equation to properly specify the problem.

My earlier presumption was that the term dWt in the equation implicitly accounts for it (I read somewhere in the documentation) which seems to be true in which case the modification of including sqrt(dt) is incorrect.

Can I use for loop to check the solution of the DiffEq the way you have done in other cases?

10

Some more testing result.
I tried the solution with sigma = 0.09 which satisfies 2alphabeta >=sigma^2 but the ensemble seem to break due to negative rates. The for loop does not throw any negatives.

11

I meant that X is non negative but it is positive under the above stated assumption. There are tons of papers discribing how to simulate X on the web…

12

It definitely does.

α=0.2;σ=0.09;u₀=0.04;β=0.041
f(u,p,t) = α*(β-u)
g(u,p,t) = σ*sqrt(u)
dt = 1//250 #1//2^(5) # in fractional form
tspan = (0.0, 10.0)
u = 0.04
for j in 1:10000
  global u
  u = 0.04
  for i in 1:10000
    u = u + dt*f(u,nothing,nothing) + g(u,nothing,nothing)*sqrt(dt)*randn()
    if u < 0
      @show i,u
    end
  end
end

That is then probably related to the chosen dt.

13

Now that you have some guarantee of positivity, I would use g(u,p,t) = σ*sqrt(abs(u)) and add

sol = solve(prob,SOSRI(),isoutofdomain=(u,p,t)->any(x->x<0,u))

that would make it reject any negative step and adaptively increase the accuracy to stay positive. However, if the solution would truly go negative, what it will do is keep honing in on the true zero, in which case it will exit with a warning of dt=dtmin. Give it a try and see what comes up: I’m not entirely sure if with those parameters it will be negative (double check they defined the model exactly the same).

14

After making g(u,p,t)=sigma*sqrt(abs(u)) the problem disappears, even using sigma=0.3.
However adding the “isoutdomain(…)” constraint, the ensemble breaks with a warning dt=dtmin as mentioned by you.

I have found some material on the web which talks about methods to ensure positivity of CIR process and one of them is log transformation as suggested by @rveltz. I will look into these methods.

Thanks for you help.