Can this numpy code be vectorized?
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?