Questo codice numpy può essere vettorializzato?

Question

Ho scritto la seguente funzione per produrre n realizzazioni di un processo CIR per un dato insieme di parametri:

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

Questo codice funziona, ma ora mi chiedo se sia possibile rimuovere il ciclo e vettorializzarlo completamente. Ho faticato a trovare un modo per farlo, poiché l'operazione eseguita nel ciclo è ricorsiva (i valori nella colonna successiva della matrice dipendono in modo non lineare dai valori delle colonne precedenti).

Inoltre, questo codice potrebbe essere semplificato? In particolare, sento che potrebbe esserci una complessità inutile nel modo in cui accedo all'ultima colonna dell'array utilizzando

np.reshape(x[:, -1], (-1, 1))

Answer 1

Poiché l'output di ogni passaggio dipende dal passaggio precedente, è impossibile eliminare completamente il ciclo a meno che la formula non possa essere semplificata per consentire un calcolo multi-passaggio. Tuttavia, il codice può sicuramente essere ancora migliorato.

• La chiamata np.concatenateè generalmente inefficiente poiché richiede la riallocazione della memoria e la copia dei dati. Il più a lungo possibile, si dovrebbe preallocare tutta la memoria necessaria e aggiornarla iterativamente.

• Alcuni calcoli possono essere spostati fuori dal ciclo per migliorare l'efficienza. Ad esempio, la np.random.normalchiamata può essere eseguita una sola volta. Lo stesso con alpha * delta_t.

• È possibile scegliere nomi di variabili più significativi.

• È sempre una buona pratica avere docstring e suggerimenti sul tipo nel codice.

Ecco una versione del codice migliorato:

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

Non ho molta familiarità con il significato di tutti i parametri nella formula. Ci sono due cose di cui ho dubbi nel tuo calcolo:

• La tariffa \$r\$nella formula viene inizializzato con lo stesso parametro muutilizzato nel calcolo iterativo. È previsto?
• L'incremento \$dW_t\$ha utilizzato una varianza fissa di 1. Dovrebbe essere questo il caso o dovrebbe essere un altro parametro della funzione?