Can this numpy code be vectorized?

Question

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))

Answer 1

Since the output of each step depends on the previous step, it is impossible to entirely eliminate the loop unless the formula can be simplified to allow a multi-step computation. However, the code can definitely still be improved.

• Calling np.concatenate is generally inefficient since it requires memory reallocation and data copying. As long as possible, one should preallocate all the needed memory and update it iteratively.

• Some computation can be moved out of the loop to improve efficiency. E.g., the np.random.normal call can be done only once. Same with alpha * delta_t.

• More meaningful variable names could be chosen.

• It is always a good practice to have docstrings and type hints in your code.

Here is one version of improved code:

def cir_simulations(alpha: float, mu: float, sigma: float, delta_t: float, sim_steps: int, num_sims: int):
    """
    Simulate the CIR process.

    Parameters:
    Input
    ...

    Output
    ...
    """
    output_shape = (num_sims, sim_steps + 1)
    sim_results = np.empty(output_shape)
    sigma_dW = np.random.normal(0, sigma * np.sqrt(delta_t), size=output_shape) 
    alpha_dt = alpha * delta_t
    sim_results[0, :] = r_prev = mu

    for r_t, sigma_dWt in zip(sim_results[1:], sigma_dW):
        r_t[:] = r_prev = (mu - r_prev) * alpha_dt + np.sqrt(r_prev) * sigma_dWt

    return sim_results

I'm not very familiar with the meaning of all the parameters in the formula. There are two things that I have doubts in your computation:

• The rate \$r\$ in the formula is initialized with the same parameter mu used in the iterative computation. Is it expected?
• The increment \$dW_t\$ used a fixed variance of 1. Is this supposed to be the case or should it be another function parameter?