Workshop 2 - SOLUTION, Algorithms and Data Analysis

NumPy, financial time series and portfolios

Author

Alberto Dorantes, Ph.D.

Published

September 28, 2026

Abstract

This is an INDIVIDUAL workshop. We move from loops to vectorization with NumPy arrays, and we use random numbers to simulate uncertain outcomes (Monte Carlo simulation). Then we work with real financial time series: we import monthly prices, calculate simple and continuously compounded returns, holding-period returns, annualized mean returns and volatilities, and portfolio returns and risk with matrix algebra.

1 General directions

Work in a new Google Colab notebook named W2-Algorithms-YourFirstName-YourLastName, and share it with cdorante@tec.mx (Edit privileges). Follow the same rules as in Workshop 1:

Re-type and run the code; write your own notes.
Predict before you run; verify every result of a challenge (with a formula, Excel, or an assert).
If you use an LLM, write your algorithm first and paste your prompt (in quotes).

Grading: notebook (60%) + Canvas Quiz W2 (40%), due on the same date as the workshop.

We will use the following libraries. Colab already has all of them installed:

import numpy as np                  # arrays, vectorization, random numbers, linear algebra
import pandas as pd                 # data frames and time series
import matplotlib.pyplot as plt     # plots

The import ... as ... statement loads a library and gives it a short alias. From now on, np.log() means “the log function of the NumPy library”. These aliases (np, pd, plt) are a universal convention: you will see them in almost every Python program for Finance.

2 From loops to vectorization

2.1 NumPy arrays

In Workshop 1 we stored cash flows in lists. Lists are flexible, but they do not do arithmetic element by element:

A = [3, 3, 4, 6, 8]
B = [-2, -3, -1, 3, 5]
print(A + B)      # the + operator CONCATENATES lists; it does not sum them!

[3, 3, 4, 6, 8, -2, -3, -1, 3, 5]

To sum the two streams of cash flows with lists, we need a loop:

total = []
for i in range(len(A)):
    total.append(A[i] + B[i])
print(total)

[1, 0, 3, 9, 13]

A NumPy array is a collection of numbers designed for mathematics. Operations on arrays are applied element by element, without writing a loop. This is called vectorization:

A = np.array([3, 3, 4, 6, 8])
B = np.array([-2, -3, -1, 3, 5])
print(A + B)          # element-by-element sum
print(A * 2)          # every element times 2
print(A / (1.15))     # every element divided by 1.15

[ 1  0  3  9 13]
[ 6  6  8 12 16]
[2.60869565 2.60869565 3.47826087 5.2173913  6.95652174]

Internally, NumPy still runs a loop, but it is written in a fast compiled language (C), and we do not need to write it.

2.2 Example: NPV of a project with vectorization

A project has 2 products, A and B, with the following expected free cash flows (millions of pesos) for years 1 to 5. The initial investment is 10 million, and the discount rate is 15%. What is the NPV of the project?

rate = 0.15
investment = 10
years = np.arange(1, 6)                  # array([1, 2, 3, 4, 5])
project_flows = A + B                    # vectorized sum of the two products
discount_factors = 1 / (1 + rate) ** years   # one discount factor per year
pv = project_flows * discount_factors        # element by element
npv = pv.sum() - investment
print("Discount factors:", discount_factors.round(4))
print(f"NPV = {npv:.4f} million pesos")

Discount factors: [0.8696 0.7561 0.6575 0.5718 0.4972]
NPV = 4.4512 million pesos

Read the code carefully: there is no loop, but each line operates on 5 numbers at the same time. np.arange(1, 6) works like range(1, 6) but returns an array.

2.3 Why vectorization matters: speed

Let’s compare a loop and vectorized code to calculate the mean of 5 million numbers:

import time

x = np.random.default_rng(1).normal(size=5_000_000)   # 5 million random numbers

start = time.time()
total = 0
for value in x:
    total = total + value
mean_loop = total / len(x)
time_loop = time.time() - start

start = time.time()
mean_vec = x.mean()
time_vec = time.time() - start

print(f"Loop: {time_loop:.3f} seconds; vectorized: {time_vec:.4f} seconds")
print(f"The vectorized version is about {time_loop / time_vec:,.0f} times faster")

Loop: 0.699 seconds; vectorized: 0.0040 seconds
The vectorized version is about 177 times faster

Both give the same result, but the vectorized version is hundreds of times faster. With financial datasets of millions of rows, this is the difference between seconds and hours.

2.4 Useful array functions

r = np.array([0.02, -0.01, 0.03, 0.015, -0.005])
print("sum:", r.sum(), " mean:", r.mean(), " std:", r.std(ddof=1).round(5))
print("max:", r.max(), " position of the max:", r.argmax())
print("cumulative sum:", r.cumsum())
print("cumulative product of growth factors:", np.cumprod(1 + r).round(5))
print("positive returns:", r[r > 0])       # filtering with a condition

sum: 0.05  mean: 0.01  std: 0.01696
max: 0.03  position of the max: 2
cumulative sum: [0.02  0.01  0.04  0.055 0.05 ]
cumulative product of growth factors: [1.02    1.0098  1.04009 1.0557  1.05042]
positive returns: [0.02  0.03  0.015]

ddof=1 indicates that the standard deviation is calculated with N-1 in the denominator (the sample standard deviation). r[r > 0] is a boolean filter: r > 0 creates an array of True/False, and only the True positions are kept.

2.5 Two-dimensional arrays and the axis parameter

A 2D array (a matrix) has rows and columns. Its shape tells us the dimensions. Many functions accept an axis parameter: axis=0 operates down the rows (one result per column), and axis=1 operates across the columns (one result per row):

# 3 months (rows) x 2 assets (columns) of returns
R = np.array([[0.01, 0.03],
              [0.02, -0.01],
              [-0.01, 0.02]])
print(R.shape)
print("mean of each asset (axis=0):", R.mean(axis=0))
print("mean of each month (axis=1):", R.mean(axis=1))

(3, 2)
mean of each asset (axis=0): [0.00666667 0.01333333]
mean of each month (axis=1): [0.02  0.005 0.005]

2.6 Matrix algebra: the `@` operator

The return of a portfolio is the weighted average of the returns of its assets:

R_P=\sum_{i=1}^{N}w_iR_i=w_1R_1+w_2R_2+\dots+w_NR_N

This sum of products is a dot product, and in Python it is written with the @ operator:

w = np.array([0.6, 0.4])              # weights of the 2 assets
print("Portfolio return in each month:", R @ w)

Portfolio return in each month: [0.018 0.008 0.002]

R @ w multiplies each row of R (one month) by the weights, and adds the products: it calculates the portfolio return for all months in one line. The portfolio variance also has a compact matrix formula:

\sigma_P^2=w^{T}\,\Sigma\,w

where \Sigma is the variance-covariance matrix of the asset returns:

cov = np.cov(R, rowvar=False)          # covariance matrix of the columns
print(cov)
var_p = w @ cov @ w
print("Portfolio variance:", var_p, " volatility:", np.sqrt(var_p))

[[ 0.00023333 -0.00018333]
 [-0.00018333  0.00043333]]
Portfolio variance: 6.533333333333335e-05  volatility: 0.008082903768654762

In the Hedge Funds course, formulas like W @ mu_a and w @ COV @ w are everywhere. Now you know how to read them.

3 Random numbers and Monte Carlo simulation

The future is uncertain. A Monte Carlo simulation generates thousands of possible scenarios with random numbers, and then summarizes the results: the average outcome, the worst 5% of the outcomes, the probability of reaching a goal, etc.

We create a random number generator with a seed. The seed makes the random numbers reproducible: with the same seed, everybody gets the same numbers.

rng = np.random.default_rng(seed=123)
sample = rng.normal(loc=0.01, scale=0.05, size=5)   # 5 draws from a normal distribution
print(sample)

[-0.03945607 -0.00838933  0.07439626  0.01969872  0.05601154]

loc is the mean and scale is the standard deviation. size can be a number or a shape, like (1000, 12) for a matrix of 1,000 rows and 12 columns.

Example. A stock has monthly returns with a mean of 1% and a standard deviation of 5%. If I invest $100 today, what is the distribution of my wealth after 12 months?

rng = np.random.default_rng(seed=123)
n_sims = 10_000
monthly = rng.normal(0.01, 0.05, size=(n_sims, 12))   # 10,000 scenarios x 12 months
wealth = 100 * np.prod(1 + monthly, axis=1)           # compound the 12 months of each scenario

print(f"Average final wealth: {wealth.mean():.2f}")
print(f"5th percentile: {np.percentile(wealth, 5):.2f}")
print(f"Probability of losing money: {(wealth < 100).mean():.2%}")

plt.hist(wealth, bins=60)
plt.title("Simulated wealth after 12 months (10,000 scenarios)")
plt.xlabel("Wealth"); plt.ylabel("Frequency")
plt.show()

Average final wealth: 112.60
5th percentile: 84.21
Probability of losing money: 26.99%

Notice the trick in the last print: (wealth < 100) is an array of True/False, and the mean of True/False values is the proportion of True values (True counts as 1, False as 0). This gives us a probability in one line.

4 CHALLENGE 1: Retirement plan (deterministic)

You are planning to save for retirement over the next 30 years. You will invest $6,000 a month in a stock account and $3,000 a month in a bond account. The stock account is expected to return 12% annually, compounded monthly, and the bond account will pay 5% annually, also compounded monthly. When you retire, you will combine your money into an account with a 9% annual return (compounded monthly). How much can you withdraw each month from your account, assuming a 25-year withdrawal period? At the end of the 25 years, your balance must be zero. Assume that all deposits and withdrawals are made at the end of each month.

Draw a mental map (a timeline) of the problem: what happens from month 1 to month 360, and from month 361 to month 660? You can draw it on paper and paste a picture in your notebook.
Write the inputs, the output and your algorithm (pseudocode).
Write a function fv_monthly_deposits(deposit, apr, months) that calculates the future value of the monthly deposits with a loop, and verify it with the future value of an annuity formula: FV=PMT\times\frac{(1+r)^{N}-1}{r}.
Calculate the value of each account at retirement, the total, and the monthly withdrawal, using the payment of an annuity: PMT=PV\times\frac{r}{1-(1+r)^{-N}}.

4.1 SOLUTION: Challenge 1

Mental map. The problem has two stages:

Stage 1 (months 1 to 360), saving: two accounts grow in parallel. Each month, each balance earns its monthly interest and receives a deposit.
Stage 2 (months 361 to 660), withdrawing: the total is moved to a new account at 9%. The balance at the beginning of stage 2 is the present value of 300 equal withdrawals.

Then, the future value of stage 1 is the present value of stage 2. This “one stage feeds the next” logic is the key of the problem.

Algorithm.

FV_stocks = future value of 6,000 per month, 360 months, 12%/12
FV_bonds  = future value of 3,000 per month, 360 months, 5%/12
PV_retirement = FV_stocks + FV_bonds
withdrawal = payment of an annuity with PV = PV_retirement, 300 months, 9%/12

def fv_monthly_deposits(deposit, apr, months):
    """Future value of equal deposits made at the END of each month (a loop)."""
    r = apr / 12
    balance = 0
    for month in range(months):
        balance = balance * (1 + r) + deposit    # interest on what I had, plus the new deposit
    return balance

def fv_annuity(deposit, apr, months):
    """Same calculation with the closed-form formula (to verify)."""
    r = apr / 12
    return deposit * ((1 + r) ** months - 1) / r

def annuity_payment(pv, apr, months):
    """Equal end-of-month payment that brings a balance pv to zero after `months`."""
    r = apr / 12
    return pv * r / (1 - (1 + r) ** (-months))

fv_stocks = fv_monthly_deposits(6000, 0.12, 360)
fv_bonds = fv_monthly_deposits(3000, 0.05, 360)

# Verification: the loop and the formula must agree
assert abs(fv_stocks - fv_annuity(6000, 0.12, 360)) < 0.01
assert abs(fv_bonds - fv_annuity(3000, 0.05, 360)) < 0.01

total_retirement = fv_stocks + fv_bonds
withdrawal = annuity_payment(total_retirement, 0.09, 300)

print(f"Stock account at retirement: ${fv_stocks:,.2f}")
print(f"Bond account at retirement:  ${fv_bonds:,.2f}")
print(f"Total at retirement:         ${total_retirement:,.2f}")
print(f"Monthly withdrawal:          ${withdrawal:,.2f}")

Stock account at retirement: $20,969,784.80
Bond account at retirement:  $2,496,775.91
Total at retirement:         $23,466,560.70
Monthly withdrawal:          $196,930.52

Explanation. Inside the loop, balance * (1 + r) + deposit is exactly what happens in a savings account each month: last month’s balance earns interest, and then the new deposit arrives at the end of the month. After 360 iterations we get the future value. The formula gives the same result, which is our verification.

The monthly withdrawal of $196,930.52 is much larger than the monthly savings of $9,000. This is the power of compounding over 30 years: most of the retirement money is interest, not deposits (you deposited only 9{,}000\times360=\$3{,}240{,}000).

A final check: if we withdraw that amount for 300 months, the balance should end at zero:

balance = total_retirement
for month in range(300):
    balance = balance * (1 + 0.09 / 12) - withdrawal
print(f"Balance after 300 withdrawals: {balance:.6f}")

Balance after 300 withdrawals: 0.000001

The final balance is practically zero (tiny differences are due to rounding in the computer).

5 CHALLENGE 2: Retirement plan with uncertainty (Monte Carlo)

In reality, the stock account does not earn exactly 1% every month. Assume that the monthly simple return of the stock account is normally distributed with a mean of 1% and a standard deviation of 4.5%. The bond account still earns exactly 5%/12 per month.

Create a random number generator with rng = np.random.default_rng(seed=2026), and simulate a matrix of monthly stock returns with rng.normal(0.01, 0.045, size=(10000, 360)): 10,000 scenarios (rows) of 360 months (columns). Use exactly this code so that everybody gets the same numbers.
Simulate the value of the stock account at retirement for each scenario. Hint: start with an array of 10,000 zeros (np.zeros(10000)), and loop over the 360 months (columns); in each month, update all scenarios at the same time with vectorization: balance = balance * (1 + returns[:, m]) + 6000.
Add the value of the bond account (from Challenge 1), and calculate the monthly withdrawal for each scenario (your annuity_payment function works with arrays, too!).
Report: (a) the average monthly withdrawal, (b) the median monthly withdrawal, (c) the 5th percentile of the monthly withdrawal, and (d) the probability that the monthly withdrawal is lower than the deterministic withdrawal of Challenge 1.
Show a histogram of the monthly withdrawals and explain with your own words why the probability in (d) is higher than 50%, even though the average monthly return is the same 1% as in Challenge 1.

5.1 SOLUTION: Challenge 2

rng = np.random.default_rng(seed=2026)
n_sims, n_months = 10_000, 360
stock_returns = rng.normal(0.01, 0.045, size=(n_sims, n_months))

balance = np.zeros(n_sims)                   # one balance per scenario
for m in range(n_months):                    # we loop over TIME...
    balance = balance * (1 + stock_returns[:, m]) + 6000   # ...but vectorize over SCENARIOS

total_sim = balance + fv_bonds               # the bond account has no uncertainty
withdrawal_sim = annuity_payment(total_sim, 0.09, 300)    # the function works with arrays

print(f"Deterministic withdrawal (Challenge 1): ${withdrawal:,.2f}")
print(f"Average simulated withdrawal:           ${withdrawal_sim.mean():,.2f}")
print(f"Median simulated withdrawal:            ${np.median(withdrawal_sim):,.2f}")
print(f"5th percentile:                         ${np.percentile(withdrawal_sim, 5):,.2f}")
prob_lower = (withdrawal_sim < withdrawal).mean()
print(f"Probability of a lower withdrawal:      {prob_lower:.2%}")

Deterministic withdrawal (Challenge 1): $196,930.52
Average simulated withdrawal:           $196,343.92
Median simulated withdrawal:            $160,550.46
5th percentile:                         $70,078.10
Probability of a lower withdrawal:      63.72%

plt.figure(figsize=(9, 5))
plt.hist(withdrawal_sim, bins=80)
plt.axvline(withdrawal, color="red", label="Deterministic withdrawal")
plt.axvline(np.median(withdrawal_sim), color="black", linestyle="--", label="Median of the simulation")
plt.xlabel("Monthly withdrawal (pesos)"); plt.ylabel("Frequency")
plt.title("Monte Carlo simulation of the monthly withdrawal (10,000 scenarios)")
plt.legend(); plt.show()

Explanation of the code. stock_returns[:, m] selects column m of the matrix: the return of month m in the 10,000 scenarios. The line balance = balance * (1 + stock_returns[:, m]) + 6000 updates the 10,000 balances at once. We still need a loop over time, because each month depends on the previous month; but we do not need a loop over scenarios. Also note that annuity_payment was written for a single number, but since it uses only arithmetic operations, it works element by element with an array: this is vectorization for free.

Why is the probability higher than 50%? The distribution of the final wealth is skewed to the right (see the histogram): a few very lucky scenarios produce enormous withdrawals and pull the average up, while the median (the typical scenario) is lower than the deterministic result. The reason is volatility drag: with volatile returns, compounding penalizes us. For example, +10% followed by −10% leaves us with 1.10\times0.90=0.99, a loss of 1%, even though the average return is 0%. The typical (median) growth of an investment is approximately \bar{R}-\sigma^{2}/2 per period, lower than the average return \bar R. Then, in most scenarios, we end up with less money than in the deterministic plan, even though the average monthly return is the same 1%. This is why financial planners use Monte Carlo simulations instead of a single “expected return” calculation.

6 Financial time series with pandas

6.1 Importing monthly prices

We will work with monthly adjusted closing prices from December 2014 to December 2025 for 7 instruments:

Ticker	Instrument
AAPL	Apple Inc.
JPM	JPMorgan Chase (bank)
WMT	Walmart
XOM	Exxon Mobil (energy)
TSLA	Tesla
GLD	SPDR Gold Shares (an ETF that follows the price of gold)
SPY	SPDR S&P 500 ETF (an ETF that follows the S&P 500 index: our market benchmark)

The following code downloads the file only if it is not already in your Colab session, and then imports it as a data frame:

import yfinance as yf
tickers = ["AAPL", "JPM", "WMT", "XOM", "TSLA", "GLD", "SPY"]
prices = yf.download(tickers, start="2014-12-01", end="2026-01-01",
                     interval="1mo", auto_adjust=True)["Close"]
prices.index = prices.index + pd.offsets.MonthEnd(0)   # dates at the end of each month

Two important arguments of read_csv:

index_col="Date" uses the Date column as the index of the data frame (the row labels).
parse_dates=True converts the dates from text into real dates, so pandas knows that this is a time series: we can select periods, resample, and plot with a time axis.

How was this file created?

The file was downloaded from Yahoo Finance with the yfinance library. You can create your own files with any tickers (for example, Mexican stocks like WALMEX.MX or the IPC index ^MXX) with this code:

import yfinance as yf
tickers = ["AAPL", "JPM", "WMT", "XOM", "TSLA", "GLD", "SPY"]
prices = yf.download(tickers, start="2014-12-01", end="2026-01-01",
                     interval="1mo", auto_adjust=True)["Close"]
prices.index = prices.index + pd.offsets.MonthEnd(0)   # dates at the end of each month

We use a fixed file in this workshop so that everybody gets exactly the same numbers for the Canvas quiz.

Silent errors with price data

Always use adjusted prices. An adjusted price includes dividends and stock splits. With unadjusted prices, a 4-for-1 stock split looks like a −75% return!
Check the order of the columns. yfinance returns the columns in alphabetical order, not in the order you requested. Always check prices.columns before creating a vector of weights by hand.

6.2 Exploring a time series

print(prices.shape)             # rows (months) and columns (tickers)
print(prices.columns.tolist())  # names of the columns
print(prices.index[:3])         # the first 3 dates of the index

(133, 7)
['AAPL', 'GLD', 'JPM', 'SPY', 'TSLA', 'WMT', 'XOM']
DatetimeIndex(['2014-12-31', '2015-01-31', '2015-02-28'], dtype='datetime64[ns]', name='Date', freq=None)

We select data by label with .loc[rows, columns] and by position with .iloc[rows, columns]:

print(prices.loc["2020-01-31", "AAPL"])               # one value, by labels
print(prices.loc["2020-01-31":"2020-04-30", ["AAPL", "SPY"]])   # a range of dates, 2 columns
print(prices.iloc[0])                                  # the first row, by position
print(prices["TSLA"].iloc[-1])                         # the last price of TSLA

74.47563934326172
Ticker           AAPL         SPY
Date                             
2020-01-31  74.475639  292.536682
2020-02-29  65.933212  269.377838
2020-03-31  61.333584  235.740250
2020-04-30  70.863205  265.675354
Ticker
AAPL     24.403900
GLD     113.580002
JPM      45.610714
SPY     169.358154
TSLA     14.827333
WMT      23.066938
XOM      56.545490
Name: 2014-12-31 00:00:00, dtype: float64
449.7200012207031

Since prices are not comparable across stocks (a stock is not “cheaper” because its price is lower), we plot a growth index: how much $1 invested at the beginning would be worth over time. We divide each column by its first value:

growth = prices / prices.iloc[0]      # vectorized: every column divided by its first value
growth.plot(figsize=(10, 5), logy=True, title="Value of $1 invested in Dec 2014 (log scale)")
plt.ylabel("Value of $1")
plt.grid(alpha=0.3)
plt.show()

We use a logarithmic scale on the y axis (logy=True) because otherwise Tesla’s growth would make all the other lines look flat. On a log scale, equal vertical distances mean equal percentage changes.

6.3 Simple returns and continuously compounded returns

A simple return is the percentage change of the price:

R_t=\frac{P_t}{P_{t-1}}-1

A continuously compounded (cc) return, also called a log return, is the natural logarithm of the growth factor:

r_t=\ln\left(\frac{P_t}{P_{t-1}}\right)=\ln(P_t)-\ln(P_{t-1})

They are related by R_t=e^{r_t}-1 and r_t=\ln(1+R_t). For small returns they are almost equal.

R = prices.pct_change().dropna()           # simple returns
r = np.log(prices).diff().dropna()         # cc returns: difference of log prices
print(R.head(3).round(4))
print(r.head(3).round(4))

Ticker        AAPL     GLD     JPM     SPY    TSLA     WMT     XOM
Date                                                              
2015-01-31  0.0614  0.0869 -0.1254 -0.0296 -0.0846 -0.0105 -0.0544
2015-02-28  0.1008 -0.0591  0.1269  0.0562 -0.0013 -0.0124  0.0204
2015-03-31 -0.0314 -0.0215 -0.0114 -0.0157 -0.0717 -0.0141 -0.0400
Ticker        AAPL     GLD     JPM     SPY    TSLA     WMT     XOM
Date                                                              
2015-01-31  0.0596  0.0833 -0.1340 -0.0301 -0.0884 -0.0105 -0.0559
2015-02-28  0.0960 -0.0609  0.1195  0.0547 -0.0013 -0.0124  0.0202
2015-03-31 -0.0319 -0.0218 -0.0115 -0.0158 -0.0744 -0.0142 -0.0408

pct_change() calculates P_t/P_{t-1}-1 for each column.
np.log(prices) takes the log of every price, and .diff() subtracts the previous value: \ln(P_t)-\ln(P_{t-1}).
The first month has no previous price, so its return is NaN (Not a Number, a missing value). .dropna() removes that row.

Internally, pct_change() and diff() use the shift function, which moves a column one row down. You can do the same calculation explicitly:

r_aapl = np.log(prices["AAPL"] / prices["AAPL"].shift(1))
print(r_aapl.head(3))

Date
2014-12-31         NaN
2015-01-31    0.059611
2015-02-28    0.096016
Name: AAPL, dtype: float64

Why two types of returns? Each one has a property that the other does not have:

cc returns add up over time. The cc return of a year is the sum of its 12 monthly cc returns. Then, the mean, the volatility and the regressions are easier with cc returns.
simple returns add up across assets. The return of a portfolio is the weighted average of the simple returns of its assets. Then, portfolios are calculated with simple returns.

This is the same rule used in the Hedge Funds course. Learn it now!

6.4 Holding period return and annualized statistics

The holding period return (HPR) is the total return of the whole period. We can calculate it with prices, or with the sum of cc returns:

HPR=\frac{P_T}{P_0}-1=e^{\sum_{t=1}^{T}r_t}-1

hpr_prices = prices.iloc[-1] / prices.iloc[0] - 1
hpr_cc = np.exp(r.sum()) - 1
print(pd.DataFrame({"HPR (prices)": hpr_prices, "HPR (cc returns)": hpr_cc}).round(4))

        HPR (prices)  HPR (cc returns)
Ticker                                
AAPL         10.1098           10.1098
GLD           2.4893            2.4893
JPM           5.9311            5.9311
SPY           2.9953            2.9953
TSLA         29.3305           29.3305
WMT           3.7993            3.7993
XOM           1.0863            1.0863

Both columns are equal, which confirms that cc returns add up over time.

To compare assets, we annualize the mean and the volatility of monthly cc returns:

Annual cc mean return =12\times\bar r_{monthly}, and the annual expected simple return is E[R]=e^{12\bar r}-1
Annual volatility =\sigma_{monthly}\times\sqrt{12}

The volatility grows with the square root of time because the variance of a sum of independent returns is the sum of their variances: Var(12\text{ months})=12\times\sigma^2_{monthly}.

stats = pd.DataFrame({
    "Annual E[R]": np.exp(12 * r.mean()) - 1,
    "Annual volatility": r.std() * np.sqrt(12),
})
stats["Return / risk"] = stats["Annual E[R]"] / stats["Annual volatility"]
stats.round(4)

	Annual E[R]	Annual volatility	Return / risk
Ticker
AAPL	0.2447	0.2673	0.9154
GLD	0.1203	0.1420	0.8473
JPM	0.1924	0.2404	0.8005
SPY	0.1342	0.1495	0.8974
TSLA	0.3637	0.5709	0.6371
WMT	0.1533	0.1871	0.8190
XOM	0.0691	0.2618	0.2641

Notice that r.mean() and r.std() calculate one number per column: pandas applies the function to each column automatically, as NumPy does with axis=0.

7 CHALLENGE 3: Risk and return of individual assets

Using the monthly cc returns r:

Create a function summary_stats(returns) that receives a data frame of monthly cc returns and returns a data frame with one row per asset and the following columns: HPR, annual expected return (e^{12\bar r}-1), annual volatility, the return/risk ratio, the worst month (minimum monthly return), and the percentage of months with a negative return.
Sort the table by annual volatility (method .sort_values()). Which asset had the highest HPR? Which one had the highest return/risk ratio? Are they the same asset?
Explain with your own words what information the HPR hides, using the example of TSLA.
Using .loc, calculate the HPR of each asset only for the year 2020 (from the price of December 2019 to the price of December 2020). Which asset had the best 2020?

7.1 SOLUTION: Challenge 3

def summary_stats(returns):
    """Summary statistics of monthly cc returns: one row per asset."""
    table = pd.DataFrame({
        "HPR": np.exp(returns.sum()) - 1,
        "Annual E[R]": np.exp(12 * returns.mean()) - 1,
        "Annual volatility": returns.std() * np.sqrt(12),
        "Worst month": returns.min(),
        "% negative months": (returns < 0).mean(),
    })
    table["Return / risk"] = table["Annual E[R]"] / table["Annual volatility"]
    return table

stats_table = summary_stats(r).sort_values("Annual volatility")
stats_table.round(4)

	HPR	Annual E[R]	Annual volatility	Worst month	% negative months	Return / risk
Ticker
GLD	2.4893	0.1203	0.1420	-0.0873	0.4773	0.8473
SPY	2.9953	0.1342	0.1495	-0.1334	0.3030	0.8974
WMT	3.7993	0.1533	0.1871	-0.1698	0.4091	0.8190
JPM	5.9311	0.1924	0.2404	-0.2544	0.4167	0.8005
XOM	1.0863	0.0691	0.2618	-0.3036	0.4697	0.2641
AAPL	10.1098	0.2447	0.2673	-0.1999	0.4242	0.9154
TSLA	29.3305	0.3637	0.5709	-0.4578	0.4621	0.6371

Explanation. The function builds a data frame from a dictionary: each key is a column name, and each value is a pandas Series with one number per asset (the tickers become the index). For the percentage of negative months, (returns < 0) produces a data frame of True/False, and .mean() gives the proportion of True values in each column, as in the Monte Carlo example.

best_hpr = stats_table["HPR"].idxmax()
best_ratio = stats_table["Return / risk"].idxmax()
print(f"Highest HPR: {best_hpr}; highest return/risk ratio: {best_ratio}")

Highest HPR: TSLA; highest return/risk ratio: AAPL

idxmax() returns the label (here, the ticker) of the maximum value, not the value itself.

What does the HPR hide? The HPR only compares the first and the last price. It tells nothing about the path: how volatile the investment was, how deep the losses were in between, or how many months were negative. TSLA has by far the highest volatility of the group: an investor in TSLA had to live through several months with losses of more than 20%, and drops of more than 50% from its peaks. An investor who needed the money in one of those bad moments would have realized a very different return. This is why we always analyze return and risk together.

HPR in 2020

hpr_2020 = prices.loc["2020-12-31"] / prices.loc["2019-12-31"] - 1
print(hpr_2020.sort_values(ascending=False).round(4))

Ticker
TSLA    7.4344
AAPL    0.8231
GLD     0.2481
WMT     0.2332
SPY     0.1833
JPM    -0.0553
XOM    -0.3621
dtype: float64

prices.loc["2020-12-31"] selects the row of that date (all tickers). Dividing two rows is vectorized: one HPR per ticker. We could also sum the cc returns of 2020: np.exp(r.loc["2020"].sum()) - 1, where r.loc["2020"] selects all the months of 2020 (pandas understands partial dates in a time-series index).

8 Portfolios with matrix algebra

A portfolio is a set of weights that add up to 1. With the matrix of monthly simple returns R (one column per asset), the monthly portfolio returns are R @ w:

ASSETS = ["AAPL", "JPM", "WMT", "XOM", "GLD"]     # we leave out SPY (the benchmark) and TSLA
w_eq = np.repeat(1 / len(ASSETS), len(ASSETS))    # equally weighted: 20% each
print(w_eq, w_eq.sum())

port_R = R[ASSETS] @ w_eq          # monthly simple returns of the portfolio
port_r = np.log(1 + port_R)        # convert them to cc returns to calculate statistics
print(port_r.head(3))

[0.2 0.2 0.2 0.2 0.2] 1.0
Date
2015-01-31   -0.008437
2015-02-28    0.034728
2015-03-31   -0.023971
dtype: float64

Selecting the columns with R[ASSETS] guarantees that the order of the columns matches the order of the weights, which is exactly the silent error we warned about above.

9 CHALLENGE 4: Diversification and random portfolios

Calculate the annual expected return and annual volatility of the equally weighted portfolio port_r (use the same formulas as in Challenge 3).
Compare the volatility of the portfolio with the weighted average of the volatilities of the 5 assets. Which one is smaller? Explain with your own words why (hint: look at the correlation matrix r[ASSETS].corr()).
Simulate 5,000 random long-only portfolios of the 5 assets: generate the weights with rng = np.random.default_rng(seed=42), W = rng.random((5000, 5)), and normalize each row so that it adds up to 1: W = W / W.sum(axis=1, keepdims=True). For each portfolio calculate its annual expected return and annual volatility with matrix algebra:
- mu = np.exp(12 * r[ASSETS].mean()).values - 1 (vector of annual expected returns)
- COV = (12 * r[ASSETS].cov()).values (annual covariance matrix)
- expected returns of all portfolios: W @ mu
- volatility of portfolio k: np.sqrt(W[k] @ COV @ W[k]) (use a loop or a list comprehension)
Plot the portfolios in a risk-return scatter plot and mark the individual assets. Report the weights of the portfolio with the highest return/risk ratio.

9.1 SOLUTION: Challenge 4

eq_ret = np.exp(12 * port_r.mean()) - 1
eq_vol = port_r.std() * np.sqrt(12)
avg_vol = w_eq @ (r[ASSETS].std() * np.sqrt(12))

print(f"Equally weighted portfolio: E[R] = {eq_ret:.4f}, volatility = {eq_vol:.4f}")
print(f"Weighted average of the individual volatilities = {avg_vol:.4f}")
print(r[ASSETS].corr().round(2))

Equally weighted portfolio: E[R] = 0.1738, volatility = 0.1305
Weighted average of the individual volatilities = 0.2197
Ticker  AAPL   JPM   WMT   XOM   GLD
Ticker                              
AAPL    1.00  0.29  0.27  0.19  0.07
JPM     0.29  1.00  0.20  0.52 -0.16
WMT     0.27  0.20  1.00  0.08  0.11
XOM     0.19  0.52  0.08  1.00 -0.08
GLD     0.07 -0.16  0.11 -0.08  1.00

The volatility of the portfolio is much lower than the weighted average of the volatilities of its assets. This is diversification: the assets do not move perfectly together (their correlations are much lower than 1, and gold is almost uncorrelated with the stocks), so the bad months of some assets are partially compensated by the good months of others. The portfolio variance formula w^T\Sigma w captures this through the covariances. Only if all correlations were exactly 1 would the portfolio volatility equal the weighted average of the volatilities.

Random portfolios

rng = np.random.default_rng(seed=42)
W = rng.random((5000, len(ASSETS)))
W = W / W.sum(axis=1, keepdims=True)      # each row adds up to 1

mu = np.exp(12 * r[ASSETS].mean()).values - 1
COV = (12 * r[ASSETS].cov()).values

port_ret = W @ mu                                              # 5,000 expected returns at once
port_vol = np.array([np.sqrt(w @ COV @ w) for w in W])         # list comprehension over the rows
ratio = port_ret / port_vol

best = ratio.argmax()
print("Best return/risk ratio:", ratio[best].round(4))
print(pd.Series(W[best], index=ASSETS).round(4))

Best return/risk ratio: 1.4445
AAPL    0.2090
JPM     0.1723
WMT     0.1566
XOM     0.0112
GLD     0.4508
dtype: float64

Explanation.

rng.random((5000, 5)) creates uniform random numbers between 0 and 1. Dividing each row by its sum (axis=1) makes the weights add up to 1. keepdims=True keeps the sums as a column (shape 5000 × 1), so that the division is applied row by row (this is called broadcasting).
W @ mu multiplies a (5000 × 5) matrix by a vector of 5 expected returns: we get the 5,000 portfolio expected returns in one line.
For the volatility, we use a list comprehension: for w in W iterates over the rows of W (one portfolio at a time), and np.sqrt(w @ COV @ w) is the formula \sqrt{w^T\Sigma w}.
ratio.argmax() returns the position of the portfolio with the highest ratio, which we use to get its weights.

plt.figure(figsize=(9, 6))
sc = plt.scatter(port_vol, port_ret, c=ratio, cmap="viridis", s=6)
plt.colorbar(sc, label="Return / risk")
asset_vol = np.sqrt(np.diag(COV))
plt.scatter(asset_vol, mu, color="red", marker="D")
for i, a in enumerate(ASSETS):
    plt.annotate(a, (asset_vol[i], mu[i]), xytext=(6, 0), textcoords="offset points")
plt.scatter(port_vol[best], port_ret[best], color="black", marker="*", s=200, label="Best return/risk")
plt.xlabel("Annual volatility"); plt.ylabel("Annual expected return")
plt.title("5,000 random long-only portfolios")
plt.legend(); plt.grid(alpha=0.3); plt.show()

Most individual assets (red diamonds) lie inside or below the cloud of portfolios: for almost any single asset, there is a portfolio with more return for the same risk, or less risk for the same return. The upper-left border of the cloud is the efficient frontier, which you will study formally in the Portfolio and Hedge Funds courses. np.diag(COV) extracts the variances (the diagonal of the covariance matrix), and enumerate gives us the position i and the ticker a to label each asset.

10 CHALLENGE 5: Code reading

Predict the output before running each piece of code. Then run it and explain any difference.

5.1

a = np.array([1, 2, 3, 4])
print(a * a, (a > 2).sum(), a[1:3].mean())

5.2

x = np.array([[1, 2], [3, 4], [5, 6]])
print(x.shape, x.sum(axis=0), x.sum(axis=1))

5.3 The following code should calculate the annual volatility from monthly returns. What is wrong?

annual_vol = r["AAPL"].std() * 12

5.4

p = pd.Series([100, 110, 99], index=["Jan", "Feb", "Mar"])
print(p.pct_change().round(2).tolist())

10.1 SOLUTION: Challenge 5

5.1 a * a is element by element: [1 4 9 16]. (a > 2) is [False False True True], and its sum counts the True values: 2. a[1:3] is [2, 3], whose mean is 2.5.

5.2 x has 3 rows and 2 columns: shape (3, 2). sum(axis=0) adds down the rows (one result per column): [9 12]. sum(axis=1) adds across the columns (one result per row): [3 7 11].

5.3 The volatility grows with the square root of time. It must be r["AAPL"].std() * np.sqrt(12). Multiplying by 12 overstates the annual volatility by a factor of 12/\sqrt{12}\approx3.46. (The mean is multiplied by 12; the standard deviation by \sqrt{12}.)

5.4 The first month has no previous value: nan. Then 110/100-1=0.10 and 99/110-1=-0.10: [nan, 0.1, -0.1]. Notice that the price went up 10% and then down 10%, but it did not return to 100: this is the volatility drag of Challenge 2.

a = np.array([1, 2, 3, 4])
print(a * a, (a > 2).sum(), a[1:3].mean())
x = np.array([[1, 2], [3, 4], [5, 6]])
print(x.shape, x.sum(axis=0), x.sum(axis=1))
p = pd.Series([100, 110, 99], index=["Jan", "Feb", "Mar"])
print(p.pct_change().round(2).tolist())

[ 1  4  9 16] 2 2.5
(3, 2) [ 9 12] [ 3  7 11]
[nan, 0.1, -0.1]

11 DataCamp online course

You MUST TAKE Chapter 3: Arrays in Python from the course Introduction to Python for Finance.

12 W2 submission

Submit the link of your Google Colab notebook in Canvas, and make sure that you shared it with cdorante@tec.mx (Edit privileges).
Answer the Canvas Quiz W2 before the deadline, with your notebook open.