# 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?

You probably have to create separate GPU arrays for $V_{t+1}$ and $V_t$.