ito process as numerical solution of SDE

[yizhang-yiz]

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

Current Version:

v2.19.1

[bob-carpenter]

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]

@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]

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]

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

[syclik]

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]

@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]

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]

@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]

@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.

[yizhang-yiz]

I'm closing this as indicated in discourse discussion(https://discourse.mc-stan.org/t/ito-process-as-numerical-solution-of-stochastic-differential-equation/9192/36?u=yizhang) the way Stan samples doesn't support the scheme to be implemented using the proposed function. Further investigation required to support SDE inference through a different approach(State space model, for example).