jlperla / knowledgediffusionsimulations.jl Goto Github PK

• Post a separate issue with your code that shows SSAStepper() drops the diffusive part (and it works fine without it). Independent of our code.

Post an example called "Simulating a large number of independent jumps" to the JumpDiffEq repo.

• Get the baseline code with the AffectIndex approach, maximum, etc. working
• Same with using a closures, and jump.... Benchmark both to see if there any speed difference.
  Impression was that with JumpProblem(prob, DirectFW()), jump...) etc. it wouldn't matter. Speed difference on N = 300.
• Post the better version on the DiffEqJump repository
• Figure out the ensemble printing See https://docs.juliadiffeq.org/latest/tutorials/sde_example.html#Ensemble-Simulations-1
  • Plot the mean and the max.

• Add trajectories to the params
• Add in the jump_algorithm = SRIW1() to the parameters
• Move rate and μ and σ and affect! into the params. Still refer back to p for the parameter choice.
  • Should be able to rename them without the _SDE.
• Single simulation: generate a new p with N = 4 and plot them all, change the title.
• Do we want a function to generate the jump problem given a p
  • create_jump_problem(p)
    Basically the stuff in

x_iv = rand(p.iv_dist, p.N)  # draws the inital condition
prob = SDEProblem(μ_SDE, σ_SDE, x_iv ,(0.0, p.t[end]), p)
rate(u,p,t) = p.β*p.N
jump = ConstantRateJump(rate,affect!)
jump_prob = JumpProblem(prob,Direct(),jump)

• Make a create_ensemble_problem(p) which calls create_jump_problem
• Try to strip down DifferentialEquations into DiffEqJump, StochasticDiffEq, DiffEqCallbacks and a few others.
• Flip over to EnsembleDistributed() and see if it works (keeping a single process)

Wewant to simulate a system of stochastic ODE where we try to limit the sizes of stochastic steps to ensure that it doesn't growth faster than a particular speed. This is a little ad-hoc, and that there is no simple markovian stochastic process that this would converge to, but I wonder if anyone has ever done anything along these lines? If the variable is in log-space, then to ensure that the growth rate is not more than g, for a given dt we would want to make sure that the shocks don't go larger than g dt. With that, a discretized version of http://docs.juliadiffeq.org/latest/features/noise_process.html#Direct-Construction-Example-1 could be

W(dt) ~ min( N(0, dt), g dt)

Which I think I can implement as something like

g = 0.05
function WHITE_NOISE_DIST_TRUNCATED(g, W, dt)
  if typeof(W.dW) <: AbstractArray
    return min.(g * dt, sqrt(abs(dt))*wiener_randn(size(W.dW)))
  else
    return min.(g * dt, sqrt(abs(dt))*wiener_randn(typeof(W.dW))
  end
end
TruncatedWienerProcess(g, t0,W0,Z0=nothing) = NoiseProcess(t0,W0,Z0, (W, dt) -> WHITE_NOISE_DIST_TRUNCATED(g, W, dt);kwargs)

We should try this out as the alternative in the http://docs.juliadiffeq.org/latest/tutorials/sde_example.html#Example-4:-Colored-Noise-1 interface.

@PooyaFa We decided to replace the VariableRateJump with a ConstantRateJump approximation due to a bug in the underlying code for that. See the DiffEq docs:

For example, if there are multiple jumps whose rates only change when one of them occur, than that set of jumps is a constant rate jump. If a jump's rate depends on the differential equation, time, or by some value which changes outside of any constant rate jump, then it is denoted as variable.

And:

(Note: if your rate is only "slightly" dependent on the solution of the differential equation, then it may be okay to use a ConstantRateJump. Accuracy loss will be related to the percentage that the rate changes over the jump intervals.)

So basically, we're now assuming the jump rates are constant in between jumps (i.e., we refresh all the jump rates whenever any firm gets a jump). Instead of before, where we updated the jump rates continuously over time (i.e., with respect to u).

The way we did this is by changing the types (VariableRateJump => ConstantRateJump), and updating the constructor in the JumpSet.

The next step is to answer the following questions:

• For a given simulation, can we quantify percent change in jump rates in between the jumps? I think one solution is: simulate the thing without any jumps, and use a callback to calculate what the jump rates would be at each callback point (say tracking for one particle, in a 3 particle setup.) We can then plot that and see if it's sufficiently small, bigger around the beginning/end, etc.

• If we increase the jump multiplier (rho_max), the approximation should be more accurate, because there are more jumps overall.

• If we decrease the absolute mean of the drift process/the volatility coefficient, the approximation should be more accurate.

• Can we have "phantom jumps" (maybe using a coupled jump problem?) that don't do anything, but let us adjust the approximation accuracy?

I played around with (4) and it looks possible... anyways, let me know if there are any issues. I think (1)-(3) are the most important sanity checks.

For many, many equations we really need to use a different interface where a time-step is not introduced for every jump. The natural candidate is the RegularJump in https://docs.juliadiffeq.org/latest/types/jump_types.html#Defining-a-Regular-Jump-1 which uses tau leaping.

What this does is fix the rate for a given time interval, and then count the number of jumps that occur in that interval (which is a poisson process). It then applies the number of jumps using a Stoichiometry matrix.

The problem is that this approach doesn't work so well for our process with the max(...) jumps because it has a fixed-size jump per process.

Maybe there is a way to have the modify the regular jump problem to support this?

We want to solve the SDE with the following features:

• Start from an arbitrary initial condition of N particles
• μ is the drift of the Brownian motion
• σ is the variance
• β is the arrival rate of jumps
• Upon an arrival of a jump, the go to max(z, z') where z' is drawn from the distribution.

To simulate, we will "solve" a jump-diffusion SDE. See http://docs.juliadiffeq.org/latest/tutorials/jump_diffusion.html

In the simplest version, we can solve this with constants for all of the parameters. Later we can have these a function of the state to simulate a "poor man's lucas-moll".

To get this working, see http://docs.juliadiffeq.org/latest/tutorials/jump_diffusion.html#Jump-Diffusion-1 for more details and http://docs.juliadiffeq.org/latest/tutorials/sde_example.html#Example-2:-Systems-of-SDEs-with-Diagonal-Noise-1 for solving the SDE on its own.

But it would be something like the following to setup the SDE (without the jumps):

using DifferentialEquations, Distributions,  StochasticDiffEq, DiffEqJump

function μ_SDE(du,u,p,t)
  du .= p.μ
end

function σ_SDE(du,u,p,t)
  du .= p.σ
end

p = (μ = 0.01, σ = 0.1, β = 0.1, N = 5) # if all constant
T = 10.0  # maximum time length
dist_iv = Normal() # initial condition distribution, though could hardcode
x_iv = rand(dist_iv, p.N)  # just draws from the inital condition

prob = SDEProblem(μ_SDE, σ_SDE, x_iv ,(0.0, T), p)

Now, to turn this into a jump problem it will take a little more work figuring out the interface. The issue is that we want N independent jump procesess. So maybe something like:

function affect!(integrator, n)
    integrator.u[n] = max(integrator.u[n], integrator.u[rand(1:integrator.p.N)])  # choose random integrator form 1 to N
end
jumps = [ConstantRateJump(0.2, (integrator) -> affect!(integrator, n) for n in 1:p.N)]  # add one for each n
jump_prob = JumpProblem(prob, Direct(), jumps)   # might be wrong...
sol = solve(jump_prob)
plot(sol) # Can't plot if too many points, instead use ensemble.

The thing I am unclear of here is the jumps setup, since we need N of them. Should that be an array? Should it be a JumpSet? Don't know, but ask for help on slack if having issues.

After you get the basics working, you should be able to crank up N, use http://docs.juliadiffeq.org/latest/features/ensemble.html if you want to get fancy, etc. to get fancy, etc.

Bug posted here: SciML/DifferentialEquations.jl#485

Might require asking for a little help on fixing a recipe?

This is a check you can add to https://github.com/jlperla/KnowledgeDiffusionSimulations.jl/tree/master/test, but we don't need it as a notebook for people to explore.

Take the single jump process with a fixed rate in the https://github.com/jlperla/KnowledgeDiffusionSimulations.jl/blob/master/notebooks/ensemble.jmd and then construct the equivalent thing with many jumps of the same fixed rate (using https://github.com/jlperla/KnowledgeDiffusionSimulations.jl/blob/master/notebooks/ensemble_variable_jump_rate.jmd as the base)

Start from the same initial condition, with a reasonably large number of firms (N = 100?)

With this, calcualate the statistics of the growth rates/etc. to make sure that these two things see to be roughly equivalent. For that, I would run a large number of simulations for each, get out the mean (of the means) and compare. Growth rates might be a little bit more volatile, but should also be compared.

Instead of an IID draw, we want the ability to distort the CDF you draw with.

To do this:

• Turn the current u into an empirical distribution and gets its CDF, call it F
• Transform the CDF by making it F^kappa for some parameter kappa > 0
• Draw from that F^kappa distribution.

How to do this? Basically I think you need to create the CDF of iid draws for the N particles, take its power, and then draw a uniform number to see where it lands.

Something like this, which you should test carefully.

N = 10
κ = 0.5

x = rand(N)
sorted_x = sort(x)

F = range(0.0, (N-1)/N, length = N) 
draw_quantile = rand()  # draw quantile on (0,1)
sorted_x[findlast(q -> q < draw_quantile, F.^κ)]

A good way to test code like this is to try with N=2 and make sure that things line up for simple cases.

For this distortion, implement it as an affect for both the multiple jumps and single jump process examples.

Once we get #6 working, we can try to speed this up by saving the sorted_x and F.^κ at the dt saveat points and draw from that. It should speed things up significantly to now have to calculate that each time, even if the approximation can be crude.

When we do that, we only want to keep track of the last sorted values and F.^κ and replace it with each dt. Otherwise it would take up a huge amount of space.

A few things on that notebook to fix up:

• I would reorganize to better fit into the template for the main notebook, https://github.com/jlperla/KnowledgeDiffusionSimulations.jl/blob/master/notebooks/ensemble.jmd for consistency (e.g. the functions, etc.) where you just put in what is necessary for the new noise process.
  • Start with that template and copy things in instead of the other way around.
• Set mu = 0.0 as the baseline. No need for drift.
• Rename the current g to be gamma (i.e. γ) throughout the notebook, in parameters, etc. See the next one as well:
• I think we need to tweak the shock process to undo the \sigma because otherwise our γ gets rescaled by it
  • Try something like this

function INPLACE_WHITE_NOISE_DIST_Trunc(p, rand_vec,W,dt,rng)
  DiffEqNoiseProcess.wiener_randn!(rng,rand_vec)
  rand_vec .*= sqrt(abs(dt))
    rand_vec .= clamp.(rand_vec, -p.g / p.σ  * abs(dt),  p.g/ p.σ * abs(dt)) 
end

• Note that I changed the interface to pass in p instead of the g or γ directly.

• After you get the above changes to the γ working, make sure the baseline value is large enough so that it dampens things a bit, but doesn't become quite so constant as it currently is. I suspect that the change above will help a lot, though. The choice is arbitrary for this.

• I think you can remove the equivalent to https://github.com/jlperla/KnowledgeDiffusionSimulations.jl/blob/master/notebooks/truncated_brownian_motion.jmd#L29-L67 and just let the users play around with the parameters.

The fear right now is that for some reason the p.moments data is being trampled or mixed from multiple ensembles.

A simple check is to modify the save_func (once, just as a check) to make sure that integrator.p.moments starts out with a length of 0... i.e. it doesn't start with data in it.

Creates an ensemble where we use the saveat=t and ignore our aggregation callback.
Same thing as the single-simulation, but with ensemble.
Something like

N = 4
sim = solve(ensemble_prob, SRIW1(), EnsembleSerial(), trajectories = 10)
plot(sim)    #change the `alpha`?

Spoke with Chris Rackauckas:

tl;dr

if it is negative (read: if we can confirm it's not jumping correctly, which is what we suspect anyway), I wouldn't expect a fix to this any time soon. We'd need to dedicate a big amount of time to this

Chris Rackauckas:shrug: 8:07 AM
variable rate jumps are just a very hard problem mathematically
it transforms the problem for each jump into an integral
makes a new term in the SDE
with a certain random initial condition
and then jumps whenever it hits zero
so there's like a few thousand phantom terms here trying to all exactly do a ContinuousCallback

Chris Rackauckas:shrug: 8:08 AM
and if one of them ever fails to trigger and goes negative, you might get this behavior
which has a probability of happening due to floating point error
so :shrug: it might just be really really hard

Chris Rackauckas:shrug: 8:09 AM
in the solution there is sol.u[i].u and sol.u[i].u_extended or something like that
check the second one for negative values.

From discussions with Isaac, it sounds like we are more likely to scale up with the approach using closures rather than a custom type - which is easier to read anyways. Make sure to use DirectFW() instead of Direct() there.

To do this, remove the custom callables and structs, and convert over to closures created with a comprehension.

But first do a quick benchmark on a solution to compare them with N = 200 or something like that. If you can't get that high, then go smaller as required. For that benchmarking, make sure to just do it on a single solution rather than an entire ensemble.

Instead of drawing iv = rand(N) or whatever, we can make it more flexible.

Instead, replace that says x_iv = rand(p.N) with

α = 2.0  # since in logs, shape of exponential distribution (i.e. cdf(x) = 1 - e^(-α x))
iv_dist = Exponential(1/α)  # use other distributions, or make deterministic
x_iv = rand(iv_dist, p.N)

Trying out an NaN for the first one, or something like that... We want to drop that g(0) in some way.

• Make sure the g divied by step(p.t)
• Add in an example for displaying a plot of a single simulation, not creating an ensemble
  • In the solve use saveat=p.t and you don't need the callback.
  • Can use use a plot(sol) there, maybe try it with a alpha to make it prettier.
  • Try this with an "interact" around it.
• Set up the code in a way that Interact with @manipulate works. For the single simulation and ensemble separately.
  • Use the variable naemes in the simulate.
• Try to get the code as compact as possible
• Convert to a .jmd

Jesse will post on the initial condition, but there will be an additional parameter or two.

Looking at #17, it seems that the issues with the huge negative shocks are present also in the single case (i.e., no aggregation.)

So the issue is with how we are setting up the variable jump problem itself. Something must have been lost in translation to the new closure structure.

I am assigning myself and @PooyaFa (since you wrote this code, it would be a huge help if you could help debug. But no worries if you are busy.)

After #6 we should have the data returned from the solution as only the summary statistics evaluated as some grid of spacing dt (which is not going to be every timestep used in the solution.

With this, we want to use the Ensemble features to plot out statistics form a large number of simulations of the setup.

See https://docs.juliadiffeq.org/latest/tutorials/sde_example.html#Ensemble-Simulations-1 and https://docs.juliadiffeq.org/latest/features/ensemble.html and https://docs.juliadiffeq.org/latest/features/ensemble.html#EnsembleSummary-1 for some documents.

Keep in mind that what we want to plot is not the underlying solution to the equations, but rather summary statistics (over the ensemble of simulations) of the summary statistics (over the actual "particle" solutions)! (e.g. what is the max of the min's, the 75th quantile of the means, etc.)

• For this, the first thing to do is triple check that the data being returned from the solution to the diffeq has only the summary statistics and is not returning the whole solution for every particle.
• Then try to use the out-of-the-box calculations of the ensemble and its plotting recipes to see if it correctly knows how to handle the "saveddata" stuff from the callbacks rather than the underlying solutions to the ODEs. This might be automatic, or it might be difficult.
  • There are Plots.jl recipes for the ensemble stuff, so use them whenever possible.
• Figure out how EnsembleSummary works. We will want to use that in the plots instead of displaying all simulations. See https://docs.juliadiffeq.org/latest/tutorials/sde_example.html#Ensemble-Simulations-1
• You probably want to go using Distributed; addprocs(10); or something like that first in order to use the default EnsembleDistributed() parallelization.
• What we really want to get for plots are various moments of the mean, min, g, max, and max/min of the distribution. Use whatever stuff you can get out of the box for the ensemble and ensemblesummary stuff. Maybe just showing the quantiles as in the https://docs.juliadiffeq.org/latest/tutorials/sde_example.html#Ensemble-Simulations-1 example is good enough...

• Get NBGitViewer links on the front page,
• Installation instructions (copy Perla1)
• remove binder

Start with N particles each with arrival rate of rho then have a single jump process with an arrival rate of N * rho.
The affect for this, randomly draws one of the particles.

function affect!(integrator)
  N = length(integrator.u)
  n = rand(1:N)
  n2 = rand(1:N)
  u[n] = max(u[n], u[n2])
end

Try to solve this for a larger N. (e.g. thousands) and look at the ensemble average.

See #8

The code should be based on the single jump process approach rather than having separate jumps. i.e. use https://github.com/jlperla/KnowledgeDiffusionSimulations.jl/blob/master/notebooks/ensemble.jmd as the base.

• add median to moments
• print stuff in g when growth is crazy (positive or negative)
• double-check results are sane with ConstantRateJump type
• increase T
• finish off notebook with different permutations (and closures, DirectFW), trying to pin down craziness

See https://github.com/jlperla/KnowledgeDiffusionSimulations.jl/blob/master/src/variable_rate_jumps.jl

Assigning this to myself so I don't forget to do it.

Make sure we really trust #16

In a separate notebook.

The idea is:

• Use a single solution, not an ensemble
• Set T to be VERY large
• Burn the first XXX of them
• Plot a histogram of g and display a few statistical properties.

Potentially play with a small number of agents (e.g. N=3) and run for two different initial conditions. The key test is whether the distribution is the same.

Additionally,

• Run an ensemble and then look at the T = 400 or whatever growth rate distribution across the ensemble, and plot the histogram/etc.

We want to allow for arrival rates for equation equation to be different, and a function of the u(t) and/or the integrator's last calculations. I think we can try this with something like what we did for the Jump Index. Here is a sketch using http://docs.juliadiffeq.org/latest/tutorials/jump_diffusion.html#Variable-Rate-Jumps-1

struct RateIndex{F1, F2}
  rate_index::F1
  index::F2
end
function (p::RateIndex)(u, p, t)
    return p.affect_index(u, p, t, p.index)
end

Then to use it, something like

rate_index(u, p, t, index) = 2.0  # although we would later use the index
affect_index!(integrator, index) = (integrator.u[index] = max(integrator.u[index], integrator.u[rand(1:integrator.p.N)]) )
jumps = [VariableRateJump(RateIndex(rate_index, i)), AffectIndex(affect_index!, i)) for i in 1:p.N]
jump_prob JumpProblem(prob,Direct(),JumpSet((), jumps, nothing, nothing))

As for the rate_index function to test with. We could start with the something that is a fixed number at the minimum of support, and 0 at the top. e.g.

p.ρ_max = 2.0
function rate_index(u, p, t, index)
   u_max = maximum(u)
   u_min = minimum(u)
   return p.ρ_max* (u[index] - u_max)^2 / (u_max - u_min)
end

We don't really want to save all of the points at all times. All we need is to save a reduction (and certainly not at all time periods).

To do this, see http://docs.juliadiffeq.org/latest/features/callback_library.html#SavingCallback-1

What we want to do is provide a callback to this which calculates the following,

using StatsBase
function solution_summary(u, t, integrator)
    summ = summarystats(u)

    # should hav a special case that if at first time-step, then make `g = 0` or something like that.
    prevmean = integrator.GETLASTSAVEDVALUEMEAN
    prevt= integrator.GETLASTSAVEDT

    # calculate the growth rate by looking at the saved values from the last time
    g = (summ.mean - prevmean)/(t - prevt) # already in logs, so just need the slow.
    return (g = g, mean = summ.mean, max = summ.max, min = summ.min, median = sum.median, q25 = summ.q25, q75 = summ.q75) # or whatever... faster ways to do all at once.
end

# Note, could figure out the return type of the saved values by artificially constructing a named tuple.
saved_values = SavedValues(Float64, Tuple{Float64,Float64,Float64 }) # or a named tuple?
saveat_range = 0.0:0.01:10.0
cb = SavingCallback(solution_summary, saved_values, saveat = saveat_range)

After that, we need to make sure that the Ensemble functions work on the savedvalues as opposed to the actual values. We also want to make sure that it isn't saving the actual u anywhere...

See https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/ensemble/basic_ensemble_solve.jl#L121

and others. Seems like the batch_func is used for both the Serial and Distributed... both of which are typical ways to use julia. Maybe the threads should just be ignored since the implementation is preliminary.

Worst case, use prob_func to do a copy of the data yourself.