I've written the following function to produce n realizations of a CIR process for a given set of parameters:
def cir_simulations(alpha, mu, sigma, delta_t, n, num_sims): x = np.reshape(np.array([mu] * num_sims), (-1, 1)) for i in range(0, n): x = np.concatenate((x, np.reshape(x[:, -1], (-1, 1)) + alpha * ( mu - np.reshape(x[:, -1], (-1, 1))) * delta_t + sigma * np.sqrt( np.reshape(x[:, -1], (-1, 1))) * np.sqrt(delta_t) * np.random.normal(0, 1, size=(num_sims, 1))), axis=1) return x
This code works, but I am now wondering whether it would be possible to remove the loop and fully vectorize it. I've struggled to find a way of doing this, as the operation being performed in the loop is recursive (the values in the next column of the matrix are non-linearly dependent on the previous columns values).
Additionally, could this code be simplified? In particular, I feel like there may be needless complexity in the way that I am accessing the last column of the array using
np.reshape(x[:, -1], (-1, 1))