Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ito process as numerical solution of SDE #1263

Open
yizhang-yiz opened this issue Jun 12, 2019 · 9 comments

Comments

Projects
None yet
3 participants
@yizhang-yiz
Copy link
Contributor

commented Jun 12, 2019

Description

ito process as numerical solution of stochastic differential equation, solved by potentially various of schemes. The easiest ones are Euler-Maruyama and its variants.

Example

The function generating ito process could be like this

  /*
   * Ito process by Euler scheme @c ito_process_tamed_euler_stepper, assuming system size n and m wiener processes are involved in diffusion, and there are q iid steps of wiener process.
   *
   * @tparam F1 functor type for autonomous drift function mu, with signature Eigen::Vector=F1(const Eigen::Vector& y, const std::vector& theta), both return and y are of size n.
   * @tparam F2 functor type for autonomous diffusion sigma, with signature @c Eigen::Matrix=F2(const Eigen::Vector& y, const std::vector& theta), y is of size n, return is n x m.
   * @tparam T0 current solution type, @c Eigen::Vector of size n
   * @tparam Tw type for iid R.Vs with standard normal
   *            distribution that serve each step of the Wiener
   *            process is based upon, @c Eigen::Matrix of size m x q.
   * @tparam T1 type of the parameters @c theta1 that @c F1 depends on.
   * @tparam T2 type of the parameters @c theta2 that @c F2 depends on.
   * @param f1 drift function
   * @param f2 diffusion function, its return must be left-multiplicable
   *           to @c w
   * @param y0 initial condition.
   * @param w Matrix of iid R.Vs for Wiener process, with
   *          size m x q, and the number of cols(q)
   *          defines the number of points of the process to be returned.
   * @param theta1 drift function parameters
   * @param theta2 diffusion function parameters
   * @param t final time to be reached through @c stepper.
   * @return solution matrix with size n x q.
   *         Solution i, i=0...q-1 is at time @c t/q*(i+1).
   */
  template<typename F1, typename F2, typename T0, typename Tw, typename T1, typename T2>
  inline Eigen::Matrix<typename stan::return_type<T0, Tw, T1, T2>::type, -1, -1>  
  ito_process_euler(const F1& f1, const F2& f2,
                    const Eigen::Matrix<T0, -1, 1>& y0,
                    const Eigen::Matrix<Tw, -1, -1>& w,
                    const std::vector<T1>& theta1,
                    const std::vector<T2>& theta2,
double t) 

Expected Output

An n x q matrix that contains the simulated ito process, with n being the system size(see above) and q being the time steps for the process. This matrix can be param or data, depending on the input types of initial condition y0, Wiener process step normals w, and functor parameters theta1 & theta2.

An example is the stochastic SIR model https://www.sciencedirect.com/science/article/pii/S0378437105001093, where the solution gives an Ito process of each compartment
c7f887db604853df9a2a4dd27512eda64c0a5c66

Current Version:

v2.19.1

@bob-carpenter

This comment has been minimized.

Copy link
Contributor

commented Jun 12, 2019

Thanks. @syclik---can you look this over and see if it's something we want to add? I'm afraid I don't know anything about stochastic diff eqs.

@yizhang-yiz

This comment has been minimized.

Copy link
Contributor Author

commented Jun 12, 2019

@syclik for some background there are two examples from https://discourse.mc-stan.org/t/ito-process-as-numerical-solution-of-stochastic-differential-equation/9192
demo the potential application. I'm working on the 3rd example.

@yizhang-yiz

This comment has been minimized.

Copy link
Contributor Author

commented Jun 12, 2019

The implementation used to run the examples in the above thread can be found at
https://github.com/yizhang-yiz/math/blob/feature_sde_tamed_euler/stan/math/rev/mat/functor/ito_process_integrator.hpp

Three examples are included in the unit test with fixed Brownian motion path: geometric Brownian motion, stochastic SIR model, stochastic volatility model(Heston model):
https://github.com/yizhang-yiz/math/blob/feature_sde_tamed_euler/test/unit/math/rev/mat/functor/ito_process_test.cpp

@syclik

This comment has been minimized.

Copy link
Member

commented Jun 14, 2019

Thanks. Looks like the Discourse thread is picking up steam. I'll comment there.

@syclik

This comment has been minimized.

Copy link
Member

commented Jun 14, 2019

Actually, some things that could be clarified on this issue itself. In the function signature:

  • what is the function signature of F1 and F2?
  • what's the dimension of y0, weiner_normal, theta1, and theta2?
  • why are theta1 and theta2 declared as real[] instead of vector?
  • can every argument, excluding the functors, be both parameters and data? (stan::math::var and double)

Is the primary use of this construct for convenience or for performance?

My take: if we can get the design straightened out, I think it's fine adding it as a function. I think that's the highest priority and we should focus on that.

On the implementation side, we need to make sure it's maintainable and there are enough tests to catch errors. and to help with future debugging. If there are dependencies to other libraries, we'll need to determine if it's appropriate if it requires additional dependencies. But I don't see why we shouldn't include something like this if it's designed well enough; it doesn't have to be perfect and hold forever. But it shouldn't just be slapped together either.

@yizhang-yiz

This comment has been minimized.

Copy link
Contributor Author

commented Jun 14, 2019

@syclik I've updated signature to c++, as I was posing in stan signature.

Is the primary use of this construct for convenience or for performance?

Convenience.

@yizhang-yiz

This comment has been minimized.

Copy link
Contributor Author

commented Jun 14, 2019

Allow me to redact the previous answer. The current proposed scheme(Euler) is rather easy to do in Stan, which is why the simple answer is "convenience". On the other hand, Euler is of order 1/2 for strong approximation, and naturally we may want something with higher order and/or more stable(like what happened in ODE solvers) in the future, say Milstein schemes, therefore laying out first steps is more than convenience.

@syclik

This comment has been minimized.

Copy link
Member

commented Jun 17, 2019

@yizhang-yiz: thanks for updating the top level comment. I don't think there's anything conceptually stopping adding a function like this. From a design perspective, at the issue level, I think we should walk through the function signature checking for:

  • whether the function signature looks natural: order of arguments
  • names of arguments
  • doc
  • expectation on the functors: do they throw? do they check for good arguments at all? is it on the integration function to figure it out?

And then at the implementation level, it'll be whether there's a sensible software architecture in place (if it needs to do anything tricky) and whether it can be maintained (just good coding practice + tests so we can move it later if need be).

@yizhang-yiz

This comment has been minimized.

Copy link
Contributor Author

commented Jun 21, 2019

@syclik:

whether the function signature looks natural: order of arguments

I'll leave it to you guys as it involves personal taste.

names of arguments

I was following ODE signature, except w(a matrix of iid standard normals).

doc

Maybe citing the literature for the scheme?

expectation on the functors: do they throw? do they check for good arguments at all? is it on the integration function to figure it out?

Since F1 and F2 are user-defined, I think it depends on the user. But we do want to check their returns to make sure consistent size.

From implementation level, the design is to break the functionality into stepper and ito_process. stepper maps to a specific integration scheme(Euler in this case, but likely more in the future) that outputs result at step i+1 given that in step i. ito_process uses stepper and given initial condition and wiener process to calculate the entire trajectory. One may be tempted to adopt boost odeint framework but IMO for the simple explicit scheme here it's an overkill.

I've also added some general background for SDE on discourse, shooting for common denominator of the audience.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.