import numpy as np # arrays, vectorization, random numbers, linear algebra
import pandas as pd # data frames and time series
import matplotlib.pyplot as plt # plotsWorkshop 2, Algorithms and Data Analysis
NumPy, financial time series and portfolios
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. 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).
When you finish the Workshop, download it as Jupyter Notebook (as a .ipynb file), and then submit that .ipynb file through Canvas.
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:
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 pesosRead 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.951 seconds; vectorized: 0.0059 seconds
The vectorized version is about 160 times fasterBoth 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 conditionsum: 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.T @ cov @ w
# w.T is the Transpose of the w vertical vector, so w.T it is the weight horizontal vector
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.008082903768654762In portfolio management, 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}}.
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 withrng.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_paymentfunction 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.
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:
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
WarningSilent 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.
yfinancereturns the columns in alphabetical order, not in the order you requested. Always checkprices.columnsbefore 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 TSLA74.47566986083984
Ticker AAPL SPY
Date
2020-01-31 74.475670 292.536713
2020-02-29 65.933205 269.377930
2020-03-31 61.333599 235.740234
2020-04-30 70.863228 265.675354
Ticker
AAPL 24.403896
GLD 113.580002
JPM 45.838268
SPY 169.358154
TSLA 14.827333
WMT 23.066938
XOM 56.545490
Name: 2014-12-31 00:00:00, dtype: float64
449.7200012207031Since 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.0408pct_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.059612
2015-02-28 0.096016
Name: AAPL, dtype: float64Why 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.
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.0863Both 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?
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: float64Selecting 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.
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() * 125.4
p = pd.Series([100, 110, 99], index=["Jan", "Feb", "Mar"])
print(p.pct_change().round(2).tolist())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.