Kann dieser Numpy-Code vektorisiert werden?

Question

Ich habe die folgende Funktion geschrieben, um n Realisierungen eines CIR-Prozesses für einen bestimmten Parametersatz zu erstellen :

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

Dieser Code funktioniert, aber ich frage mich jetzt, ob es möglich wäre, die Schleife zu entfernen und vollständig zu vektorisieren. Ich habe mich bemüht, einen Weg zu finden, da die in der Schleife ausgeführte Operation rekursiv ist (die Werte in der nächsten Spalte der Matrix hängen nicht linear von den Werten der vorherigen Spalten ab).

Könnte dieser Code zusätzlich vereinfacht werden? Insbesondere habe ich das Gefühl, dass die Art und Weise, wie ich mit auf die letzte Spalte des Arrays zugreife, unnötig komplex sein kann

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

Answer 1

Da die Ausgabe jedes Schritts vom vorherigen Schritt abhängt, ist es unmöglich, die Schleife vollständig zu beseitigen, es sei denn, die Formel kann vereinfacht werden, um eine mehrstufige Berechnung zu ermöglichen. Der Code kann jedoch definitiv noch verbessert werden.

• Das Aufrufen np.concatenateist im Allgemeinen ineffizient, da es eine Neuzuweisung des Speichers und das Kopieren von Daten erfordert. So lange wie möglich sollte der gesamte benötigte Speicher vorab zugewiesen und iterativ aktualisiert werden.

• Einige Berechnungen können aus der Schleife verschoben werden, um die Effizienz zu verbessern. Beispielsweise kann der np.random.normalAnruf nur einmal getätigt werden. Gleiches gilt für alpha * delta_t.

• Es könnten aussagekräftigere Variablennamen gewählt werden.

• Es ist immer eine gute Praxis, Dokumentzeichenfolgen und Tipphinweise in Ihrem Code zu haben.

Hier ist eine Version des verbesserten Codes:

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

Ich bin nicht sehr vertraut mit der Bedeutung aller Parameter in der Formel. Es gibt zwei Dinge, an denen ich Zweifel habe:

• Die Rate \$r\$in der Formel wird mit demselben Parameter initialisiert, der mubei der iterativen Berechnung verwendet wird. Wird es erwartet?
• Das Inkrement \$dW_t\$verwendete eine feste Varianz von 1. Soll dies der Fall sein oder sollte es ein anderer Funktionsparameter sein?