Last active
January 13, 2019 22:59
-
-
Save kailaix/c68b358a7b3f21e06ef084d1e23a2a15 to your computer and use it in GitHub Desktop.
CIR Model
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| using LinearAlgebra | |
| using PyPlot | |
| using SparseArrays | |
| # https://uu.diva-portal.org/smash/get/diva2:1148496/FULLTEXT01.pdf page 8 | |
| # https://regulation.pstat.ucsb.edu/sites/sa-wcmf6/5-Peng_ECIR_Model_Qidi.pdf page 13 | |
| # params | |
| r0 = 0.5 | |
| L = 4.0 | |
| a = 1.0 | |
| b = 1.0 | |
| σ = 1.0 | |
| @show σ^2 - 2*a*b | |
| N = 100 | |
| r = LinRange(0, L, N+1)[2:end] | |
| Δr = r[2]-r[1] | |
| T = 0.3 | |
| G = 0.01 | |
| NT = 100 | |
| Δt = T/NT | |
| t = LinRange(0, T, NT+1) | |
| # assembling | |
| function add(i, j, v) | |
| if i<1 || i>N || j<1 || j>N | |
| return | |
| end | |
| push!(ii, i) | |
| push!(jj, j) | |
| push!(vv, v) | |
| end | |
| function add_derive(i, c, order) | |
| if order==1 | |
| add(i, i+1, c/(2Δr)) | |
| add(i, i-1, -c/(2Δr)) | |
| elseif order==2 | |
| add(i, i+1, c/Δr^2) | |
| add(i, i-1, c/Δr^2) | |
| add(i, i, -2c/Δr^2) | |
| end | |
| end | |
| function assemble() | |
| return sparse(ii,jj,vv,N,N) | |
| end | |
| function clear!() | |
| global ii, jj, vv | |
| ii = Int[] | |
| jj = Int[] | |
| vv = Float64[] | |
| end | |
| clear!() | |
| for i = 1:N-1 | |
| add(i, i, 1/Δt) | |
| add_derive(i, (a*(b-r[i])-σ^2)/2, 1) | |
| add_derive(i, (-σ^2*r[i]/2)/2, 2) | |
| add(i, i, -a/2) | |
| end | |
| add(N, N, σ^2-a*(b-r[end])) | |
| add(N, N, σ^2*r[end]/2/Δr) | |
| add(N, N-1, -σ^2*r[end]/2/Δr) | |
| A = assemble() | |
| clear!() | |
| for i = 1:N-1 | |
| add(i, i, 1/Δt) | |
| add_derive(i, -(a*(b-r[i])-σ^2)/2, 1) | |
| add_derive(i, -(-σ^2*r[i]/2)/2, 2) | |
| add(i, i, a/2) | |
| end | |
| B = assemble() | |
| U = zeros(N, NT+1) | |
| U[:,1] = 1/sqrt(2π*G^2)*exp.(-(r .- r0).^2/2G^2) | |
| U[:,1] = U[:,1]/sum(U[:,1])/Δr | |
| for i = 1:NT | |
| U[:,i+1] = A\(B*U[:,i]) | |
| @show sum(U[:,i+1])*Δr | |
| end | |
| mesh(t, r, U) |
Author
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I implemented a simple CIR model using assembler programming.