CIR Simulation on GPUs: Arrayfun Limits and CUDA Kernels
Summary
The document asks how to simulate the Cox–Ingersoll–Ross (CIR) variance process efficiently on a GPU in MATLAB. It presents a nonnegative Euler discretization, where each step uses the prior variance, mean reversion, and a Brownian shock, then clips negative values to zero. Its sample implementation loops over time and stores paths across Monte Carlo simulations.
The reported issue is that MATLAB’s GPU arrayfun does not support the indexed, recursive updates used by the loop, while calling the function directly runs but is slower than the CPU in the author’s setup. Replies suggest that arrayfun has restrictions and that CUDA or PTX kernels may be appropriate for more flexible GPU execution. Another suggestion is to represent current and next time-step values in separate GPU arrays and consult MATLAB examples. The discussion gives implementation guidance rather than benchmark evidence or a tested solution; performance will depend on the hardware, MATLAB version, and implementation.
Key ideas
- The Euler scheme updates each CIR variance value recursively from the previous time step.
- Clipping the Euler update at zero keeps simulated variance nonnegative.
- The indexed loop cannot be directly expressed through the described GPU arrayfun call.
- CUDA kernels or separate arrays for adjacent time steps are suggested workarounds, not demonstrated solutions.
Tags
Full text
# How to simulate a CIR process using GPU and Matlab
# How to simulate a CIR process using GPU and Matlab
I am trying to simulate a CIR process using Matlab and my GPU for effeciency. At the moment i run into some implementation problems due to the recursive nature of the discretization.
The sheme I currently use is the simple Euler
$$V_t \approx \max \{V_{t-1} + \kappa_V(\bar{V}-V_{t-1})\Delta h + \sigma_V\sqrt{V_{t-1}}\sqrt{\Delta h}\mathcal{N}(0,1),0\}. $$
To proper simulate on the GPU I need some help function, like discribed in this Matlab guide
```
function [ V ] = simulation_fun( V, W, N, h, kappa_V, V_bar, sigma_V)
t=1;
while t < N % a path has N steps
V(:,t+1)=max(V(:,t)+kappa_V*(V_bar-V(:,t))*h + sigma_V*sqrt(V(:,t))*sqrt(h)*W(:,t),0);
t=t+1;
end
end
```
Than I define the model parameters
```
N=1000;
M=1000; % Number of MC simulations
h=0.01;
kappa_V=1;
V_bar=1;
sigma_V=0.2;
V=gpuArray(ones(M,N));
W=gpuArray.randn(M,N);
```
After that I call the simulation via
```
V = arrayfun(@simulation_fun,V, W, N, h, kappa_V, V_bar, sigma_V);
```
and get the follwing error
```
Error using gpuArray/arrayfun
Indexing is not supported. error at line: 3 .
```
However, if I use
```
V = simulation_fun(V, W, N, h, kappa_V, V_bar, sigma_V);
```
it works, but the simulation is roughly 2.5 times slower than only using the CPU.
Does anyone know how to correctly implement the simulation?
## Answer by ALFRAM (score 2)
https://quant.stackexchange.com/a/16780
You should write some kernel functions in CUDA (Nvidia language) for your matlab code. Arrayfun is quite restrictive and not appropriate. Look at this link http://fr.mathworks.com/help/distcomp/run-cuda-or-ptx-code-on-gpu.html for more details about matlab and parallel computing.
## Answer by jaamor (score 1)
https://quant.stackexchange.com/a/16776
There are some restrictions to using arrayfun. You can read the restrictions here.
Judging from the error, you cannot use indexes the way you are.
You probably have to create separate GPU arrays for $ V_{t+1} $ and $V_t$.
I suggest that you find similar examples in Matlab's website and try to replicate its functionality. Here is an article with information on how to program the second order wave equation on a GPU with Matlab.Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.