# Simulating CIR model with Differential Equation Package

**URL:** <https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569>\
**Category:** General Usage\
**Tags:** question\
**Created:** [September 9, 2019, 5:28pm UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569 "2019-09-09T17:28:50Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ajaychat3](https://avatars.discourse-cdn.com/v4/letter/a/ecd19e/32.png) [@Ajaychat3](https://discourse.julialang.org/u/Ajaychat3)\
**Post date:** [September 9, 2019, 5:28pm UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/1 "2019-09-09T17:28:50Z")

</div>

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

Please help me resolving this issue.

```julia
# 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
```

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [September 9, 2019, 5:51pm UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/2 "2019-09-09T17:51:00Z")

</div>

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)`

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 9, 2019, 6:39pm UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/3 "2019-09-09T18:39:39Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Ajaychat3](https://avatars.discourse-cdn.com/v4/letter/a/ecd19e/32.png) [@Ajaychat3](https://discourse.julialang.org/u/Ajaychat3)\
**Post date:** [September 10, 2019, 2:09am UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/4 "2019-09-10T02:09:15Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 10, 2019, 2:33am UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/5 "2019-09-10T02:33:47Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Ajaychat3](https://avatars.discourse-cdn.com/v4/letter/a/ecd19e/32.png) [@Ajaychat3](https://discourse.julialang.org/u/Ajaychat3)\
**Post date:** [September 10, 2019, 3:02am UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/6 "2019-09-10T03:02:25Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 10, 2019, 3:36am UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/7 "2019-09-10T03:36:06Z")

</div>

> [@Ajaychat3](#):
>
> I didn’t get you when you say stddev is much larger than 1.

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.

> [@Ajaychat3](#):
>
> If I change sqrt(u) to sqrt(u\*dt) in the equation even with sigma=0.2, it seems to work.

That’s not your original equation.

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

```julia
α=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

```

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 10, 2019, 3:43am UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/8 "2019-09-10T03:43:58Z")

</div>

This seems to be a pretty well-known result:

> **[Cox–Ingersoll–Ross model](https://en.wikipedia.org/wiki/Cox%E2%80%93Ingersoll%E2%80%93Ross_model)**
>
> In mathematical finance, the Cox–Ingersoll–Ross (CIR) model describes the evolution of interest rates. It is a type of "one factor model" (short-rate model) as it describes interest rate movements as driven by only one source of market risk. The model can be used in the valuation of interest rate derivatives. It was introduced in 1985 by John C. Cox, Jonathan E. Ingersoll and Stephen A. Ross as an extension of the Vasicek model.
> The CIR model specifies that the instantaneous interest rate 
>   
> ...

[http://janroman.dhis.org/finance/IR/CIR%20advanced.pdf](http://janroman.dhis.org/finance/IR/CIR%20advanced.pdf)

Positivity is guaranteed when

```julia
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.

---

<div class="post-metadata">

**Author:** ![Ajaychat3](https://avatars.discourse-cdn.com/v4/letter/a/ecd19e/32.png) [@Ajaychat3](https://discourse.julialang.org/u/Ajaychat3)\
**Post date:** [September 10, 2019, 4:29am UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/9 "2019-09-10T04:29:58Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Ajaychat3](https://avatars.discourse-cdn.com/v4/letter/a/ecd19e/32.png) [@Ajaychat3](https://discourse.julialang.org/u/Ajaychat3)\
**Post date:** [September 10, 2019, 5:10am UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/10 "2019-09-10T05:10:22Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [September 10, 2019, 5:40am UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/11 "2019-09-10T05:40:47Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 10, 2019, 5:57am UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/12 "2019-09-10T05:57:21Z")

</div>

> [@Ajaychat3](#):
>
> The for loop does not throw any negatives.

It definitely does.

```julia
α=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`.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 10, 2019, 6:00am UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/13 "2019-09-10T06:00:23Z")

</div>

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

```julia
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).

---

<div class="post-metadata">

**Author:** ![Ajaychat3](https://avatars.discourse-cdn.com/v4/letter/a/ecd19e/32.png) [@Ajaychat3](https://discourse.julialang.org/u/Ajaychat3)\
**Post date:** [September 10, 2019, 8:56am UTC](https://discourse.julialang.org/t/simulating-cir-model-with-differential-equation-package/28569/14 "2019-09-10T08:56:18Z")

</div>

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.