Solving difference equation

General Usage Performance
,
1

This is my code to solve a system of difference equations:

function trajectoryDiscrete(sys_eq, init_cond, NIter, par1, par2)
    neq = length(init_cond)
    traj= zeros(Float64,NIter+1,neq)
    traj[1,:] = init_cond
    sys(x) = sys_eq(x, par1, par2);
     for i in 1:NIter
        init_cond = sys(init_cond)
        traj[i+1,:] = init_cond
    end
    return traj
end

Here is a sample system of difference equations:

function test_difference_eq(x,r,k)
    a = 5.0; mu = 0.5; d = 0.2; 
    dx = x[1]*exp(r*(1-x[1]/k)-x[2]/(a+x[1]^2))
    dy = x[2]*exp(mu*x[1]/(a+x[1]^2)-d)
    return @SVector [dx, dy]
end

Here is the result:

tr = @time trajectoryDiscrete(test_difference_eq, [0.5, 0.3], 10000, 1, 2)

  0.000414 seconds (3 allocations: 156.484 KiB)
10001×2 Matrix{Float64}:
  0.5             0.3
 13.2216          0.380422
  0.000982828  2586.64
  ⋮
  0.0             2.5e-323
  0.0             2.5e-323

Query: If I remove @SVector from the return statement of test_difference_eq then it gives

0.000655 seconds (10.00 k allocations: 1.068 MiB)
10001×2 Matrix{Float64}:
 0.5      0.3
 0.99971  0.257598
 1.5792   0.229228
 ⋮
 2.0      2.5e-323
 2.0      2.5e-323

Ah! you guessed it right, I am new to Julia, trying to settle here ditching MATLAB. My queries are:

  1. What magic does @SVector do so that allocation is reduced. Can you give me an idea when, where and how I can use @SVector and @SMatrix?

  2. Can this code be optimized further?

My ultimate goal is to solve the system for a range of r=range(1, 5, length = 501) and k=range(2, 4, length = 501) in nested for loop, store last few point of tr to find its period using another user-defined function seqper (see my previous post) and finally store each r-k pair along with the associated period in a text file.

What are your suggestions?

Any form of help is appreciated.

Blockquote

2 Likes
2

SVector is allocated on stack when it’s small, think of it as a NTuple basically, it’s especially useful when you need to use a small matrix / vector over and over again.

2 Likes
3

Your code is row major. You may want to take a look at the 18.337 course notes which include optimizing a difference equation as one of its early lectures.

5 Likes
4

@MathGuy SolveDSGE.jl solves infinite horizon systems of stochastic & deterministic difference equations

@ChrisRackauckas will SciML ever be able to solve difference equations? Or is it only intended for continuous time?

1 Like
5

See DynamicalSystems.jl. It solves dicrete problmes like above. Just use trajectory. For further details see this link.

6

DynamicalSystems.jl uses the DiffEq DiscreteProblem interface for that.

https://diffeq.sciml.ai/stable/types/discrete_types/

It’s not the most talked about one because, well, it’s not so hard to write the loops, but it exists and puts everything into the same interface as the DiffEq solutions. It’s mostly used as the basis for the discrete jump problems, but you can also just directly define and solve a discrete dynamical system with it. It’s also what ModelingToolkit.jl’s DiscreteSystem targets.

2 Likes
7

@ChrisRackauckas the simplest system economists solve is an infinite horizon, 2nd order, boundary value (not initial value), non-linear system of difference equations:
\left[ \begin{array}{c} (c_{t})^{-\sigma } = \beta (\alpha (k_{t+1})^{\alpha-1} + (1-\delta)) (c_{t+1})^{-\sigma } \\ k_{t+1} = k_t^\alpha + (1-\delta)k_{t} - c_t \end{array}\right]
Boundary values:
Initial condition: k(0)=k_0 given
Terminal condition: \lim_{t\to \infty} \beta^t (c_{t})^{-\sigma } k_{t+1} =0

This is easy to solve w/ SolveDSGE.jl. (include the stochastic extensions…)
Can SciML handle this yet?
Does SciML want to eventually?

8

Yeah, I owe @jlperla some work getting state space models onto the SciML interface.

Is that implicit discrete?

2 Likes
9

No, but it would be nice at some point. The mixed initial value/boundary value problems here are a little non-standard for people outside of economics though, so I think having solvers that layer on top of the sciml stuff is probably better for now. The steady state transversality conditions are also very different from the standard sciml cases, but will hopefully get there.

My guess is that the first step is to get the simulation of state space models (and their differentiation) correct first, which is something that we will hopefully make progress on this fall.

4 Likes