今天，我们介绍如何用Python绘制资产组合理论（又称“均值-方差分析”）中的资产组合散点图，它的最外侧、上半部的点所连成的线被称为效前沿，efficient frontier。当资产数为2时，各个点可以用一条曲线连起来，但资产数为3及3个以上时，所有点就组成一个面。

全部程序位于最下，请首先确保安装好如下库：pandas、numpy、matplotlib、datetime和yfinance。安装好后，我们就分5步，介绍如何获取股票数据、计算每只股票的收益率的均值和它们的方差-协方差矩阵（variance-covariance matrix，VCM）、如何生成组合、计算组合的均值和方差，最后，画出散点图。

第一步：获取股票数据。

使用Yahoo Finance（yfinance）库，导入需要分析的几只股票，这里先要查到它们的股票代码（Ticker symbol），美国股票是以英文缩写表示的，比如MSFT就是Microsoft（微软公司）的股票代码，BAC代表Bank of America（美国银行）。这里，选择频率为月度，即yf.download中的interval参数的值为'1mo'，默认是日度。

下载后，数据包含了开盘价、收盘价、调整后的收盘价乃至交易量等变量。通常，我们选择调整后的收盘价（它调整了配股、分红等影响股价的因素，此处不细说）作为股票价格。然后，使用Pandas库中的pct_change函数计算收益率。

第二步：计算每只股票收益率的均值和它们的方差-协方差矩阵。

Pandas库中的函数mean、cov直接就能计算，mean函数默认就是对每一列求均值，不像Numpy等库中的均值函数那样还要写明对行还是对列求。但是，由于我们接下来要用到矩阵计算，因此，我们需要把Pandas计算出来的值通过to_numpy函数转换为向量或矩阵。

第三步：生成权重向量（矩阵）。

这是难点，为什么？首先在理论上，所谓的一个资产组合，其本质上就是一个权重向量，这是容易被遗忘或忽视的。另外，在实际操作上，我们要生成N个组合，N可以很大，比如10^4个，这不可能通过手动一个个地生成，还有，要随机地生成，因此，我们需要用到随机数生成函数。举例：先设N=10，我们现在考虑3只股票，因此我们要生成10×3的随机数矩阵，可是我们的目的是生成权重向量（矩阵），而这些随机数每一行加起来不等于1， 因此我们还要把每一行的随机数相加，用这个和去除每一个随机数，这就实现了新的随机数的每一行加和为1，符合权重向量的定义。

但是，这里面有一个很难处理的问题，就是权重为负的情况。我们知道，资产组合中，某些股票的权重是可以为负的。那么，我们只要把随机数的生成区间设置成由负到正不就可以了吗？这么做当然能生成随机数，但是，如果这样做，像我们上一段说的“加和后然后做除法”，其中有几行的加和，就可能因为正负相加的原因，出现很接近0的数，如此，分母很接近0，就会导致权重出现极端值，几千几万都可能，从而进一步导致下一步中组合的均值、方差出现极端值，最终导致画图失败。为此，解决的方法是，要么不考虑负权重，要么在这一步允许负权重，但是在下一步，删去极端值。

第四步：计算组合的期望收益率和波动率。

这是计算的难点。首先，我们必须知道，就一个组合而言，其期望收益率和方差（波动率的平方）的计算公式分别为：

这里和分别为第二步计算的股票收益率均值向量和方差-协方差矩阵。但是，这里编程的难点是，我们不止有一个权重向量，我们生成的是N个权重向量，也就是一个权重矩阵，比如记作。因此，需要动脑几分钟，想想怎样相乘才能一次性地得到N个组合的均值、N个组合的方差。如果实在想不出来，就用循环，根据上面两个式子计算每一个组合的期望收益率和方差，这样循环N次。

这里，我们采用点乘与矩阵乘相配合的方法。用权重矩阵点乘股票均值向量，再按列加和，就得到了N个组合的均值。计算组合方差时，先用矩阵乘法（符号为@），把“N*股票数”（比如股票数为3）维度的权重矩阵与“股票数*股票数”维度的方差-协方差矩阵相乘，得到“N*股票数”的一个矩阵，再用它点乘，再按列加和，就得到了N个组合的方差。最后，别忘开根号，因为我们用的是波动率。

第五步：绘图。

使用matplotlib.pyplot，绘制散点图。参数中，可以选择点的透明度、大小和colormap的样式。最后，用plt.show()呈现图像，括号中可根据自己电脑的情况，看看是否要加block=True，有时不加可能图像无法呈现。两种随机数生成方法绘制的散点图如下：

以上就是绘制资产组合图像的简要介绍，有些细节，由于时间和精力限制问题，难以一一说明，程序可直接复制运行，除了有些地方可能需要根据自己的Python版本和环境进行微调外，大概率一次就能成功。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime
import yfinance as yf

# Step 1. Download montly level stock data
start_date = datetime(2018, 1, 1)
end_date = datetime(2021, 1, 1)
tickers = ['AAPL', 'TSLA', 'MSFT', 'BAC', 'JPM']
data = yf.download(tickers, start=start_date, end=end_date, interval='1mo')
# use adjusted close prices
dataclose = data['Adj Close']  # data.columns --> MultiIndex, e.g., ([('Adj Close', 'AAPL'), ('Adj Close', 'MSFT'), ...]
# calculate returns
stockret = dataclose.pct_change()

# Step 2. Cacluate average returns and the variance-covariance matrix
retmean = stockret.mean().to_numpy()
retmean.shape = (1, len(tickers)) 
vcm = stockret.cov().to_numpy()

# Step 3. Generate random variables and weights
N = 20000  # 1000 above is fine, bigger than 10^5 makes plotting slow
rvmethod = 0  # Set it to be 0 or other values
if rvmethod == 0:
    rv = np.random.uniform(0.1, 5, size=[N, len(tickers)])
else:
    rv = np.random.normal(0, 5, size=[N, len(tickers)])  # This method allows for negative weights but need to exclude extreme values later

rvsum = np.sum(rv, axis=1)
rvsumrep = np.matlib.repmat(rvsum, len(tickers), 1).T  # replicate the column vector to facilitate division in the next line
weightmat = rv/rvsumrep  # important step: generate weight matrix
print(np.sum(weightmat, axis=1))  # check whether weights in each row sum up to one

# Step 4. Cacluate mean return and volatility of porfolios
portret = np.sum(weightmat*retmean, axis=1)  # 点乘时，dim(weightmat)=N*股票数, dim(retmean)=1*股票数，Python能自动理解是weightmat的每一行与retmean点乘。
portret.shape = (N, 1)
portvol = np.sqrt(np.sum((weightmat@vcm) * weightmat, axis=1))  # @是矩阵乘，*是点乘，即相同位置的元素相乘，这是计算的难点，实在不行就一个个地循环。
portvol.shape = (N, 1)

# this step is mainly for dropping extreme values
tab0 = pd.DataFrame(np.hstack((portvol, portret)), columns=['portvol', 'portret'])
if rvmethod == 0:
    tab = tab0
else:
    tab = tab0[tab0['portvol'] <= 0.15]

# Step 5. Draw scatter plot
f1 = plt.figure(1)
plt.scatter(tab['portvol'], tab['portret'], alpha=0.5, s=100, c=np.linspace(1, 10, len(tab)), cmap='viridis')
plt.xlabel("Portfolio Volatility", fontsize=14)
plt.ylabel("Portfolio Expected Return", fontsize=14)
plt.grid(linestyle=':')
plt.show(block=True)