Week 5: Simulation Methods

Version 5.0 - August 2026

Author

Affiliation

Jon Danielsson

London School of Economics

5 Week 5

Many risk questions have no closed-form answer. When a payoff depends on the path prices take, or when the quantity of interest is the far tail of a distribution rather than its centre, the practical approach is to simulate, drawing a large number of possible outcomes and looking at the distribution that results.

This week generates random numbers, uses them to compare distributions and to sample discrete outcomes in a credit-rating example, examines how much the answer moves as the sample grows, and packages a simulation into a function.

Before you start

Open the main course folder and check that the Prices.csv and Returns.csv files created in Week 2 are present. This week uses the distribution ideas from Week 3 and introduces no additional libraries.

By the end of this session

You should be able to explain what a random seed guarantees, compare normal and Student-t simulations, simulate discrete credit ratings, explain why larger simulated samples give more stable results, and write a function that simulates a random walk.

5.1 The plan for this week

Make random numbers
Compare distributions
Sample discrete outcomes for a credit risk example
Compare small and large simulated samples
Learn how to make functions

5.2 Loading data and libraries

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

Returns = pd.read_csv("Returns.csv", index_col="date", parse_dates=True)
Prices = pd.read_csv("Prices.csv", index_col="date", parse_dates=True)

Check your result

Both data frames should contain the ten securities created in Week 2 and share the same date index.

5.3 Random numbers

Random numbers are drawn from a distribution by calling a method on a numpy.random.Generator object. The Generator is created by np.random.default_rng(), and the draws that follow are method calls on that object.

rng = np.random.default_rng(888)

print(rng.normal(size=1))
print(rng.normal(loc=10, scale=0.5, size=10))
print(rng.standard_t(df=5, size=1))

[0.26303893]
[10.25128263  9.69822768 10.5805266  10.44903575  9.9186223   9.51590937
 10.25256047  9.87189609  9.15127419 10.84668435]
[0.00678864]

Every call above advances rng to the next values in its stream, so the numbers printed here will not repeat if the chunk is run again in the same session. What a seed guarantees is the sequence a given Generator produces from the start, not any particular number appearing at any particular place, so what follows describes distributional and statistical properties of the numbers rather than quoting specific values.

5.3.1 Seeds and reproducibility

To see what a seed actually guarantees, create three throwaway generators purely for this illustration, two seeded identically and one seeded differently.

demo_a = np.random.default_rng(888)
demo_b = np.random.default_rng(888)
demo_c = np.random.default_rng(666)

print(demo_a.normal(size=5))
print(demo_b.normal(size=5))
print(demo_c.normal(size=5))

[ 0.26303893  0.50256526 -0.60354463  1.1610532   0.8980715 ]
[ 0.26303893  0.50256526 -0.60354463  1.1610532   0.8980715 ]
[ 0.05252847  0.40854786 -0.19149297 -2.34163103  2.22114217]

demo_a and demo_b produce identical draws, since they were seeded with the same value, and demo_c produces a different sequence. Note what default_rng() returns. It is an explicit Generator object that we hold onto and pass to wherever randomness is needed, rather than a hidden global stream that every random call in the session quietly reads from and mutates. The older np.random.seed() and np.random.rand() interface works the second way and has been discouraged since numpy 1.17. The sections that follow draw from the single rng created above, and the last section shows the other arrangement, in which a function is given a seed and creates a generator of its own.

5.4 Comparing distributions using random numbers

To show the fat tails of the Student-t distribution compared to the normal, draw 1000 points from each using rng. A Student-t with three degrees of freedom has variance three, so it is scaled to variance one first. Otherwise its points would spread further for two reasons at once, and only one of them is the tails.

nu = 3
t_scale = np.sqrt((nu - 2) / nu)

rnd_normal = rng.standard_normal(1000)
rnd_t = rng.standard_t(df=nu, size=1000) * t_scale

fig, ax = plt.subplots()
ax.scatter(range(1000), rnd_t, color="red", s=15, label="Student-t")
ax.scatter(range(1000), rnd_normal, color="blue", s=15, label="Normal")
ax.set_title("Random points from a normal and Student-t distribution")
ax.legend(loc="lower right", frameon=False)
plt.show()

Repeat this chunk in the same session and the points move, because rng has already advanced past these draws. A clean render of the whole document reproduces them exactly, since rng is created again from the same seed and the same calls follow in the same order. Either way the pattern is the same. The Student-t points take on more extreme values than the normal points, visible as the red points reaching further from zero than the blue ones. Both sets of draws have the same variance, so this is a consequence of fat tails alone. The Student-t distribution assigns more probability to large deviations than the normal does, and with only three degrees of freedom that effect is pronounced.

5.5 Discrete random sampling in credit risk

Not every quantity of interest is continuous. Judging a risk model against history, as in backtesting, means resampling observations that have already happened, and a credit rating is one of a fixed set of labels rather than a number on a line.

Suppose a bond portfolio is to be simulated, each bond drawn according to the historical probability of each rating. The distributions used above cannot produce such draws, but rng.choice() can.

rng.choice(a, size, p) takes an array a, a size and optional probabilities p, returning an array of random draws.

# Credit ratings from best to worst
ratings = ["AAA", "AA", "A", "BBB", "BB", "B", "CCC", "D"]

# Illustrative probabilities for new bond issues
probs = [0.02, 0.05, 0.15, 0.35, 0.25, 0.12, 0.05, 0.01]

# Simulate ratings for 100 bonds in a portfolio
bond_ratings = rng.choice(ratings, size=100, p=probs)

Counting the outcomes with plain pd.Series.value_counts() would sort by frequency rather than by rating, which is a less natural order for a rating scale running from AAA to D. Marking the ratings as an ordered categorical keeps them in the sequence they were defined in:

rated = pd.Categorical(bond_ratings, categories=ratings, ordered=True)
rating_counts = pd.Series(rated).value_counts(sort=False)
rating_counts

AAA     2
AA      8
A      17
BBB    26
BB     29
B      13
CCC     4
D       1
Name: count, dtype: int64

How many bonds are investment grade (BBB or better)?

investment_grade = np.isin(bond_ratings, ["AAA", "AA", "A", "BBB"]).sum()
print("Investment grade bonds:", investment_grade, "out of", len(bond_ratings))

Investment grade bonds: 53 out of 100

Check your result

The probabilities should sum to 1 and the rating counts should sum to 100. The exact counts change with the random stream, but BBB and BB should usually be the largest categories because they have the largest probabilities.

The same counts as a bar chart:

fig, ax = plt.subplots()
ax.bar(rating_counts.index, rating_counts.values, color="lightblue")
ax.set_title("Simulated bond portfolio by credit rating")
ax.set_ylabel("Number of bonds")
ax.set_xlabel("Credit rating")
plt.show()

The same approach is used in credit risk modelling to simulate portfolio compositions. The draws here are independent of one another, so the simulation says nothing about the correlated downgrades that make a concentrated portfolio dangerous. That needs a model of dependence between issuers. As with the earlier draws, the counts are particular to this run and this generator. A different seed would give a different composition with the same broad shape, concentrated around BBB and BB, since those carry the largest probabilities.

5.6 Small and large sample properties

How closely a random sample follows the distribution it was drawn from depends on its size. Draw standard normal samples of four sizes and compare each with the normal density. The same bin edges are used in every panel, so the bars can be compared directly:

sample_sizes = [50, 100, 1000, 100000]
x = np.linspace(-3, 3, 100)
bins = np.linspace(-4, 4, 21)

fig, axes = plt.subplots(2, 2, figsize=(8, 8))
for i in range(len(sample_sizes)):
    ax = axes.flat[i]
    n = sample_sizes[i]
    draws = rng.standard_normal(n)
    ax.hist(draws, bins=bins, density=True, color="lightgrey")
    ax.plot(x, stats.norm.pdf(x), color="red", linewidth=2)
    ax.set_title(f"Sample size {n}")
    ax.set_ylim(0, 0.5)
fig.tight_layout()
plt.show()

With only 50 or 100 draws, the histogram is noticeably irregular and its bars depart from the smooth red density curve in places. By 1,000 draws the shape is recognisable, and with 100,000 the histogram tracks the normal density closely across almost its whole range. Each bar is the proportion of draws falling into one bin, which is a sample average, so the law of large numbers applies bar by bar. Every bar converges to the probability the density assigns to its bin as the sample grows. It is the reason later seminars can rely on large simulated samples and not on small ones.

5.7 Writing functions

A function packages a calculation so that it can be run again with different inputs. def names it and lists what it needs, the indented block is the calculation, and return hands back the result. Data and settings that vary from one call to the next should arrive through the arguments rather than be picked up from the surrounding document. Libraries imported at the top of the file are the exception.

The example below simulates a geometric random walk, a price path whose daily log returns are normal draws. What the function does is fixed. The security, the number of days, the volatility and the seed are supplied at each call.

def random_walk(name, n, sigma, seed=None):
    """Simulate and plot n days of a geometric random walk from a price of 100."""
    rng = np.random.default_rng(seed)
    steps = rng.normal(0, sigma, n)
    price = 100 * np.exp(np.cumsum(np.insert(steps, 0, 0)))
    fig, ax = plt.subplots()
    ax.plot(price)
    ax.set_title(f"Simulated {name} price over {n} days")
    ax.set_xlabel("Day")
    ax.set_ylabel("Price")
    plt.show()
    return price

The walk is three lines. rng.normal(0, sigma, n) draws n daily log returns, np.insert(steps, 0, 0) puts a return of zero in front of them so that the path opens at its starting value, and np.cumsum() followed by np.exp() turns the accumulated returns into prices beginning at 100.

Three arguments have to be supplied at every call, the security’s name, the number of days and the daily volatility. seed is given a default value of None, so it may be left out, and np.random.default_rng(None) takes its starting point from the operating system instead of from a number we choose. A seeded call therefore repeats exactly, and an unseeded one differs every time.

The volatility is taken from the course data rather than invented, so the simulated path moves as that security has actually moved:

walk = random_walk("JPM", 250, Returns["JPM"].std(), seed=888)

Calling it again with the same seed produces the same path:

walk_again = random_walk("JPM", 250, Returns["JPM"].std(), seed=888)

Leaving the seed out gives a different path on every call, and a different one each time the document is rendered:

walk_unseeded = random_walk("Citigroup", 500, Returns["C"].std())

5.8 Recap

5.8.1 In this seminar, we have covered:

Drawing random numbers from specified distributions with a numpy.random.Generator
Creating a Generator from a seed and holding onto it, rather than resetting a hidden global stream
Comparing distributions visually using random samples
Small and large sample properties of random numbers
Discrete sampling techniques for financial applications, using an ordered categorical to keep results in a natural order
Writing a function, with required arguments and one that has a default, to package a repeated calculation

5.8.2 Some new functions used:

np.random.default_rng() — create a Generator object from a seed
rng.normal() — generate random numbers from a normal distribution
rng.standard_normal() — generate random numbers from the standard normal distribution
rng.standard_t() — generate random numbers from a Student-t distribution
rng.choice() — random sampling from an array with optional probabilities
pd.Categorical() — mark values as belonging to an ordered set of categories
.value_counts() — create a frequency table
np.isin() — test membership in an array
ax.hist() — create histograms
ax.bar() — create bar charts
def — define a function, with a default value for an argument that may be left out
np.cumsum() — accumulate a sequence of values
np.insert() — place a value at a given position in an array
np.exp() — exponentiate, here turning accumulated returns into prices

5.9 Exercises

5.9.1 Apply the method

Change the credit-rating portfolio size from 100 to 1,000. Before running the code, predict what will happen to the percentages in each category. Compare the two bar charts and explain why the larger sample is normally closer to the stated probabilities.

Then give random_walk() a fourth argument, a daily drift, added to each log return before the path is accumulated. Simulate one security with a drift of zero and the same security with a small positive drift, using the same seed for both, and explain which features of the path the drift changes and which it leaves to the random draws.

5.9.2 Compare results

Change the seed and rerun the normal, Student-t and credit-rating simulations. Identify which numerical results change and which qualitative conclusions stay the same. Then display rating percentages rather than counts.

5.9.3 Extend the analysis

Simulate credit ratings separately for two industries using different probability vectors, then compare the combined portfolio with a single pooled simulation.
Write bootstrap_volatility(returns, n_boot, rng), use it to construct a 95% bootstrap interval for JPM volatility, and explain the sampling uncertainty it reveals.
Generate normal and Student-t returns matched to an actual security’s mean and standard deviation, then compare both with the real returns using the QQ plots introduced in Week 3.