Simulating for loop 100 times with a matrix

user_14656239 0
Viewed 228

I am currently trying to run a simulation of MBS a 100 times. I have written the for loop that I need, however I need to make this loop run 100 times. After consulting with a friend, I believe that I need to specify the "size" parameter in np.random.normal as a matrix, however my coding skills are limited and I would greatly appreciate your help with doing so. Specifically, for a sequence of correlation parameter ρ (rho) between 0 and 1, I need to simulate 100 MBSs, and report the average payoff of each tranche across the simulations Below is my code with notes.

EDIT: I appreciate the code proposed in the answers and it is indeed very useful. I have one last hurdle now to include: How do I include a payment structure that is sequential. Specifically, The junior tranche is the first to absorb losses from the underlying collateral pool and does so until the portfolio loss exceeds 5% (i.e. the proportion of defaults exceeds 10%) at which point the junior tranche becomes worthless. The mezzanine tranche begins to absorb losses once the portfolio loss exceeds 5% and continues to do so until the portfolio loss reaches 10% (i.e. the proportion of defaults exceeds 20%). Finally, the senior tranche absorbs portfolio losses in excess of 10%. As of now we've only considered an average payoff given a specific share of the total payoff.

rho_list = [0,0.2, 0.5,0.6, 1] 

# parameters 
n_borrowers = 10
payoff_default = 0.5
payoff_nodefault = 1
threshold = -1.65

# draw of income shocks (I have to draw new values of s and eps for each simulation)
s = np.random.normal(0,1,size=1) # common s for all borrowers
eps = np.random.normal(0,1,size=n_borrowers) # each borrower has their own eps

for rho in rho_list :

  # compute the borrower's income
  x = np.sqrt(rho) * s + np.sqrt(1-rho) * eps

  # which borrower defaults?
  loan_payoff = (x < threshold) * payoff_default + (x >= threshold) * payoff_nodefault

  # total pool 
  total_payoff = np.sum(loan_payoff)

  # how much does each investor receive from total_payoff?
  senior_payoff = 0.82 * total_payoff 
  mezz_payoff = 0.12 * total_payoff 
  junior_payoff = 0.6 * total_payoff 

  print('total: {}, senior: {}, mezz: {}, junior: {}'.format(total_payoff, senior_payoff, mezz_payoff, junior_payoff))

# next steps (this is what I need help with)
# repeat this for 100 simulations and compute the average payoff to each investor
# is it possible to generate the income for all simulations in one step? 
# Idea: specify the "size" parameter in np.random.normal as a matrix
finance python
1 Answers

I'm not certain I've understood the problem entirely correctly, but I think this does what you want without using any for loops by taking advantage of numpy's broadcasting. I'm by no means an expert in numpy, and multidimensional calculations are something I'm not super comfortable with but I believe my logic is sound. I'd be more than happy for any feedback.

Solution

# Setup
import pandas as pd
import numpy as np


investors = np.array([0.82, 0.12, 0.06])
rhos = np.linspace(0, 1, 11)[..., None, None]

default = 0.5
nodefault = 1
thresh = -1.65

n_sims = 100
n_borrowers = 10

s = np.random.normal(0, 1, size=(1, n_sims, 1))
eps = np.random.normal(0, 1, size=(n_sims, n_borrowers))


# Solution
x = np.sqrt(rhos) * s + np.sqrt(1 - rhos) * eps
payoffs = (x < thresh) * default + (x >= thresh) * nodefault

avgs = payoffs.sum(axis=2).mean(axis=1)
investor_payouts = avgs[..., None] * investors[None, ...]

data = np.hstack([rhos.reshape(-1, 1), investor_payouts])
df = pd.DataFrame(data, columns=["rho", "senior", "mezz", "junior"])

Output:

    rho  senior    mezz  junior
0   0.0  8.0155  1.1730  0.5865
1   0.1  7.9909  1.1694  0.5847
2   0.2  7.9991  1.1706  0.5853
3   0.3  7.9868  1.1688  0.5844
4   0.4  7.9991  1.1706  0.5853
5   0.5  7.9786  1.1676  0.5838
6   0.6  7.9745  1.1670  0.5835
7   0.7  7.9745  1.1670  0.5835
8   0.8  7.9458  1.1628  0.5814
9   0.9  7.9007  1.1562  0.5781
10  1.0  7.8310  1.1460  0.5730

Rationale

With n_sims = 100 and n_borrowers = 10, and rhos = np.linspace(0, 1, 11) we have the shapes

>>> rhos.shape
(11, 1, 1)

>>> s.shape
(1, 100, 1)

>>> eps.shape
(100, 10)

The reasoning for the shapes of rhos and s are so that broadcasting can be done more easily.

For each simulation, we essentially need to calculate the payoffs for each ρ. In essence, we want an array of shape (11, 100, 10) where along the first axis are the values of ρ, and the second and third axis are one hundred simulations of 10 borrowers.

The first term of your equation is sqrt(ρ) * s, and we want (11, 100, 1) so that we can broadcast later.

np.sqrt(rhos) * s
# shapes (11, 1, 1) * (1, 100, 1) = (11, 100, 1)

This gives us the same 100 simulated values for s, each multiplied by a different value of sqrt(ρ) (e.g., for ρ=0, which is the first value in rhos, the first row of this (11, 100) matrix is all zeros). We've added an extra dimension to get (11, 100, 1) in order to add to the second term.

The second term follows a similar logic, we want the values of sqrt(1 - ρ) to be multiplied across 100 simulations of 10 borrowers. Since eps.shape == (100, 10) and rhos.shape == (11,), and we want (11, 100, 10), we need to add two new axes to rhos:

np.sqrt(1 - rhos) * eps
# shapes (11, 1, 1) * (100, 10) = (11, 100, 10)

Now we want to combine those two terms for a final array of shape (11, 100, 10). This is why we gave the first term a new axis to get (11, 100, 1), which allows us to broadcast the values of the first term over the second term's last axis:

np.sqrt(rhos) * s + np.sqrt(1 - rhos) * eps
# shapes (11, 100, 1) + (11, 100, 10) = (11, 100, 10)

We're doing this because, in your original code, you are taking a scalar s and broadcasting it over eps, which was an array of length 10. In order to do that, numpy needed to broadcast s into an array of shape (10,) to match the shape of eps. We're doing the same thing here, except we're trying to do it for 100 simulations AND 11 different ρ values.

After all that nasty broadcasting, we arrive at an array which we can now collapse into a sum across borrowers, (total_payoff = np.sum(loan_payoff) in your original code), and then an average across all 100 simulations, which is achieved by the axis arguments to those respective functions; axis 2 has 10 elements, representing the borrowers; axis 1 has 100 elements, representing each simulation. So we use

payoffs.sum(axis=2).mean(axis=1)

Note that the calculation of the intermediary x is the same as in your original code.

At this point, we've obtained the average total payoff for 100 simulations across 10 borrowers, for 11 different values of ρ. From here we want to break out the average payoff by investor. In other words, we have 11 average payoffs (one for each ρ), and 3 investor rates, and we want to broadcast the 3 investor rates over the 11 average payoffs to get an array of shape (11, 3).

Right now avgs.shape == (11,) and investors.shape == (3,) so we need to add some axes to get our desired result:

investor_payouts = avgs[..., None] * investors[None, ...]
# shapes (11, 1) * (1, 3) = (11, 3)

Finally, the np.hstack stuff isn't necessary, that's just me stacking the ρ values with the results so that I could put everything in a dataframe. You could just as easily create the resultant dataframe in a number of other ways, depending on what you need.

user_62987120
Related