Can this numpy code be vectorized? 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))
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):
    &quot;&quot;&quot;
    Simulate the CIR process.

    Parameters:
    Input
    ...

    Output
    ...
    &quot;&quot;&quot;
    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?