Pandas-数据分析实用指南第二版-四-

Pandas 数据分析实用指南第二版（四）

原文：annas-archive.org/md5/ef72ddf5930d597094f1662f9e78e83e

译者：飞龙

协议：CC BY-NC-SA 4.0

第三部分：应用——使用 Pandas 进行真实世界分析

现在是时候将我们迄今为止所学的内容整合起来了。在本节中，我们将使用一些真实世界的数据集，进行从头到尾的分析，结合前几章所讲的所有概念，并在过程中引入一些新材料。

本节包括以下章节：

• 第七章，金融分析——比特币与股票市场

• 第八章，基于规则的异常检测

第八章：第七章：金融分析——比特币与股票市场

是时候转变思路并开始开发一个应用程序了。在本章中，我们将通过分析比特币和股票市场来探索一个金融应用。 本章将建立在我们迄今为止学到的所有知识基础上——我们将从互联网提取数据；进行一些探索性数据分析；使用pandas、seaborn和matplotlib创建可视化；计算分析金融工具性能的关键指标；并初步体验模型构建。请注意，我们这里并不是要学习金融分析，而是通过介绍如何将本书中学到的技能应用于金融分析。

本章也是本书中标准工作流程的一次突破。到目前为止，我们一直将 Python 作为一种功能性编程语言来使用。然而，Python 也支持面向对象编程，例如StockReader类（用于获取数据）、Visualizer类（用于可视化金融资产）、StockAnalyzer类（用于计算金融指标）以及StockModeler类（用于建模金融数据）。由于我们需要大量代码来使分析过程简洁且易于复现，我们将构建一个 Python 包来存放这些类。代码将在文中展示并逐步解释；然而，我们不需要自己键入或运行它——确保阅读本章的本章材料部分，正确设置环境。

本章将具有一定挑战性，可能需要多读几遍；然而，它将教会你最佳实践，并且在这里学到的技能将大大提高你的编程能力，这些技能会迅速带来回报。一个主要的收获应该是，面向对象编程（OOP）在打包分析任务方面非常有帮助。每个类应该有一个单一的目标，并且要有良好的文档。如果我们有许多类，应该将它们分布到不同的文件中，创建一个包。这样，其他人就能很容易地安装/使用它们，我们也能标准化在整个项目中执行某些任务的方式。举个例子，我们不应该让项目中的每个合作者都写自己的数据库连接函数。标准化且良好文档化的代码将节省很多后续麻烦。

本章将涵盖以下主题：

• 构建一个 Python 包

• 收集金融数据

• 进行探索性数据分析

• 对金融工具进行技术分析

• 使用历史数据建模性能

本章材料

本章中，我们将创建自己的股票分析包。这使得我们非常容易分发我们的代码，并且其他人也能使用我们的代码。该包的最终产品已上传至 GitHub：github.com/stefmolin/stock-analysis/tree/2nd_edition。Python 的包管理工具 pip 可以从 GitHub 安装包，也可以本地构建包；因此，我们可以选择以下任意方式继续操作：

• 如果我们不打算修改源代码以供个人使用，可以从 GitHub 安装。

• Fork 并克隆仓库，然后在本地机器上安装它，以便修改代码。

如果我们希望直接从 GitHub 安装，这里不需要做任何操作，因为在第一章中设置环境时已经完成了安装，数据分析入门；不过，作为参考，我们可以执行以下操作来从 GitHub 安装包：

(book_env) $ pip3 install \
git+https://github.com/stefmolin/stock-analysis.git@2nd_edition

小贴士

URL 中的@2nd_edition部分告诉pip安装标记为2nd_edition的版本。要安装特定分支上的版本，只需将其替换为@<branch_name>。例如，如果我们希望安装在名为dev的分支上开发的代码，可以使用@dev。当然，务必先检查该分支是否存在。我们还可以使用提交哈希值以相同方式抓取特定的提交。有关更多信息，请访问pip.pypa.io/en/latest/reference/pip_install/#git。

若要在可编辑模式下本地安装—即任何更改会自动在本地反映，而无需重新安装—我们使用-e标志。请在我们在第一章中创建的虚拟环境的命令行中运行以下命令。请注意，这将克隆该包的最新版本，可能与书中的版本不同（即带有2nd_edition标签的版本）：

(book_env) $ git clone \
git@github.com:stefmolin/stock-analysis.git
(book_env) $ pip3 install -r stock-analysis/requirements.txt
(book_env) $ pip3 install -e stock-analysis

重要提示

本示例使用的是通过 SSH 克隆git clone；如果尚未设置 SSH 密钥，请改为通过 HTTPS 克隆，使用如下 URL：https://github.com/stefmolin/stock-analysis.git。或者，可以按照 GitHub 上的说明首先生成 SSH 密钥。如果你有兴趣只克隆带有2nd_edition标签的版本，请参考这个 Stack Overflow 贴文：stackoverflow.com/questions/20280726/how-to-git-clone-a-specific-tag。

我们将在本章中使用此包。本书仓库中本章的目录包含我们实际分析时将使用的financial_analysis.ipynb笔记本，地址为github.com/stefmolin/Hands-On-Data-Analysis-with-Pandas-2nd-edition/tree/master/ch_07。data/文件夹包含备份文件，以防数据源自发布以来有所更改，或使用StockReader类收集数据时出现错误；如果发生这种情况，只需读取 CSV 文件并按照本章剩余内容进行操作。同样，exercises/文件夹也包含练习的备份文件。

重要说明

如果我们在使用 Jupyter Notebook 时修改了一个以可编辑模式安装的包中的文件，我们需要重新启动内核或打开一个新的 Python shell 并重新导入该包。这是因为 Python 会在导入后缓存它。其他选项包括使用importlib.reload()或 IPython 的autoreload扩展（ipython.readthedocs.io/en/stable/config/extensions/autoreload.html）。

构建一个 Python 包

构建包被视为良好的编码实践，因为它允许编写模块化代码并实现重用。使用matplotlib绘制图形时，我们不需要知道我们调用的函数内部到底在做什么——只需了解输入和输出是什么，就足以在其基础上进行构建。

包结构

window_calc.py来自第四章，聚合 Pandas DataFrame，以及viz.py来自第六章，使用 Seaborn 绘图与自定义技巧，这两个都是模块。包是由组织成目录的模块集合。包也可以被导入，但当我们导入一个包时，我们可以访问其中的某些模块，这样就不必单独导入每一个模块。这还允许我们构建可以相互导入的模块，而不必维护一个非常大的模块。

为了将模块转化为包，我们按照以下步骤进行：

1. 创建一个名为包的目录（本章使用stock_analysis）。

2. 将模块放置在上述目录中。

3. 添加一个__init__.py文件，其中包含导入包时要执行的任何 Python 代码（这个文件可以为空，也经常是空的）。

4. 在包的顶层目录（此处为stock_analysis）的同级别创建一个setup.py文件，它将向pip提供有关如何安装该包的指令。有关创建此文件的详细信息，请参见进一步阅读部分。

一旦上述步骤完成，可以使用 pip 安装该包。请注意，尽管我们的包仅包含一个目录，但我们可以根据需要构建包含多个子包的包。这些子包的创建与创建包时类似，唯一不同的是它们不需要 setup.py 文件：

1. 在主包目录（或某个其他子包内）创建一个子包目录。

2. 将子包的模块放入该目录。

3. 添加 __init__.py 文件，其中包含当导入子包时应执行的代码（此文件可以为空）。

一个包含单个子包的包的目录层次结构大致如下所示：

repo_folder
|-- <package_name>
|   |-- __init__.py
|   |-- some_module.py
|   `-- <subpackage_name>
|       |-- __init__.py
|       |-- another_module.py
|       `-- last_module.py
`-- setup.py

构建包时需要注意的其他事项包括以下内容：

• 为仓库编写README文件，以便他人了解它包含的内容（见：www.makeareadme.com/）。

• pylint 包，网址：www.pylint.org/)。

• 添加测试，确保代码修改不会破坏任何功能，并且代码能够按预期工作（请查看 pytest 包，网址：docs.pytest.org/en/latest/）。

stock_analysis 包概述

在本章中，我们将使用迄今为止讨论的各种 Python 包以及 Python 标准库，创建一个名为 stock_analysis 的 Python 包。该包位于 stock-analysis 仓库中（github.com/stefmolin/stock-analysis），其结构如下：

图 7.1 – stock-analysis 仓库的结构

我们包中的模块将包含用于进行资产技术分析的自定义类。类应为单一目的而设计，这样有助于构建、使用和调试，特别是当出现问题时。因此，我们将构建多个类，以覆盖财务分析的各个方面。我们需要为以下每个目的创建一个类：

图 7.2 – stock_analysis 包的主要主题和类

可视化包中模块之间的交互及每个类所提供的功能是很有帮助的。为此，我们可以构建统一建模语言（UML）图。

UML 图

utils.py 用于工具函数：

图 7.3 – stock_analysis 包的模块依赖关系

提示

pylint包附带了pyreverse，它可以生成 UML 图。如果已安装graphviz（www.graphviz.org/download/），从命令行运行以下命令可以生成 PNG 文件，显示模块之间的关系以及类的 UML 图（前提是已克隆代码库并安装了pylint）：pyreverse -o png stock_analysis

stock_analysis包中类的 UML 图如下所示：

图 7.4 – stock_analysis包中类的 UML 图

每个框的顶部部分包含类名；中间部分包含该类的属性；底部部分包含该类中定义的任何方法。注意从AssetGroupVisualizer和StockVisualizer类指向Visualizer类的箭头吗？这意味着这两个类都是Visualizer的一种类型。对于AssetGroupVisualizer和StockVisualizer类所显示的方法，它们在这些类中与Visualizer类中的定义不同。我们将在探索性数据分析部分更深入地探讨这一点。在本章的其余部分，我们将更详细地讨论stock_analysis包中的每个类，并利用它们的功能对金融资产进行技术分析。

收集金融数据

在第二章《使用 Pandas DataFrame》中，以及第三章《使用 Pandas 进行数据清洗》中，我们使用 API 收集数据；然而，还有其他方式可以从互联网上收集数据。我们可以使用pandas提供的pd.read_html()函数，它会为页面上找到的每个 HTML 表格返回一个 DataFrame。对于经济和金融数据，另一种选择是pandas_datareader包，它被stock_analysis包中的StockReader类用于收集金融数据。

重要说明

如果本章使用的数据源发生了变化，或者在使用StockReader类收集数据时遇到错误，可以读取data/文件夹中的 CSV 文件作为替代，以便继续跟随文本进行操作；例如：

pd.read_csv('data/bitcoin.csv', index_col='date', parse_dates=True)

StockReader 类

由于我们将在相同的日期范围内收集各种资产的数据，创建一个隐藏所有实现细节的类是有意义的，因此避免大量复制粘贴（以及潜在的错误）。为此，我们将建立StockReader类，它将使得收集比特币、股票和股票市场指数的数据变得更容易。我们可以简单地通过提供我们分析所需的日期范围创建StockReader类的实例，然后使用它提供的方法获取任何我们喜欢的数据。以下的 UML 图表提供了实现的高层次概述：

图 7.5 – StockReader 类的 UML 图表

UML 图表告诉我们StockReader类提供了一个可用股票代码（available_tickers）的属性，并且可以执行以下操作：

• 使用get_bitcoin_data()方法以所需货币拉取比特币数据。

• 使用get_forex_rates()方法拉取每日外汇汇率数据。

• 使用get_index_data()方法拉取股票市场上特定指数（如标准普尔 500 指数）的数据。

• 使用get_index_ticker()方法查找特定指数的股票市场符号（例如，在 Yahoo! Finance 上的标准普尔 500 指数的^GSPC）。

• 使用get_risk_free_rate_of_return()方法收集无风险收益率。

• 使用get_ticker_data()方法拉取股票市场上特定股票（比如 Netflix 的 NFLX）的数据。

现在我们理解了为什么需要这门课程，并对其结构有了高层次的概述，我们可以继续查看代码。由于在stock_analysis/stock_reader.py模块中有大量代码需要审查，我们将逐个部分地分解这个文件。请注意，这可能会改变缩进级别，因此请查看文件本身以获取完整版本。

模块的第一行是关于模块本身的help()，它将出现在顶部附近。这描述了我们模块的目的。接着是我们将需要的任何导入：

"""Gather select stock data."""
import datetime as dt
import re
import pandas as pd
import pandas_datareader.data as web
from .utils import label_sanitizer

注意，import语句按照PEP 8（Python 风格指南，网址为 www.python.org/dev/peps/pe…

1. 标准库导入（datetime和re）

2. 第三方库（pandas和pandas_datareader）

3. 来自stock_analysis包中另一个模块的相对导入（.utils）

在我们的导入之后，我们定义了StockReader类。首先，我们创建一个字典，将指数的股票代码映射到描述性名称_index_tickers中。注意，我们的类还有一个文档字符串，定义了它的目的。在这里，我们只会列出一些可用的股票代码：

class StockReader:
    """Class for reading financial data from websites."""
    _index_tickers = {'S&P 500': '^GSPC', 'Dow Jones': '^DJI',
                      'NASDAQ': '^IXIC'}

在构建类时，有许多特殊方法（俗称dunder 方法，因为它们的名称以双下划线开头和结尾），我们可以提供这些方法来定制类在与语言操作符一起使用时的行为：

• 初始化对象（__init__()）。

• 使对象可以进行排序比较（__eq__()、__lt__()、__gt__()等）。

• 执行对象的算术运算（__add__()、__sub__()、__mul__()等）。

• 能够使用内建的 Python 函数，例如len()（__len__()）。

• 获取对象的字符串表示，用于print()函数（__repr__()和__str__()）。

• 支持迭代和索引（__getitem__()、__iter__()和__next__()）。

幸运的是，我们不需要每次创建类时都编写所有这些功能。在大多数情况下，我们只需要__init__()方法，它会在我们创建对象时执行。（有关特殊方法的更多信息，请访问dbader.org/blog/python-dunder-methods和docs.python.org/3/reference/datamodel.html#special-method-names.)

StockReader类的对象持有数据采集的开始和结束日期，因此我们将其放入__init__()方法中。我们解析调用者传入的日期，以便允许使用任何日期分隔符；例如，我们将能够处理 Python datetime对象的输入；形如'YYYYMMDD'的字符串；或使用任何与非数字正则表达式（\D）匹配的分隔符表示日期的字符串，例如'YYYY|MM|DD'或'YYYY/MM/DD'。如果有分隔符，我们将其替换为空字符串，以便在我们的方法中使用'YYYYMMDD'格式构建日期时间。此外，如果调用者给定的开始日期等于或晚于结束日期，我们将引发ValueError：

    def __init__(self, start, end=None):
        """
        Create a `StockReader` object for reading across 
        a given date range.
        Parameters:
            - start: The first date to include, as a datetime 
              object or a string in the format 'YYYYMMDD'.
            - end: The last date to include, as a datetime
              object or string in the format 'YYYYMMDD'.
              Defaults to today if not provided.
        """
        self.start, self.end = map(
            lambda x: x.strftime('%Y%m%d')\
                if isinstance(x, dt.date)\
                else re.sub(r'\D', '', x),
            [start, end or dt.date.today()]
        )
        if self.start >= self.end:
            raise ValueError('`start` must be before `end`')

请注意，我们没有在__init__()方法中定义_index_tickers，该方法在对象创建时被调用，因为我们只需要为从这个类创建的所有对象保留一份该信息。_index_tickers类属性是私有的（按约定，前面有一个下划线），这意味着，除非该类的用户知道它的名称，否则他们不会轻易找到它（请注意，方法也可以是私有的）。这样做的目的是为了保护它（尽管不能完全保证）并且因为用户并不需要直接访问它（它是为了类的内部工作）。相反，我们将提供一个属性，我们可以像访问属性一样访问它，还会提供一个类方法，用于获取映射到字典中给定键的值。

提示

类方法是可以在类本身上使用的方法，无需事先创建类的实例。这与我们到目前为止看到的实例方法相对。实例方法是与类的实例一起使用的，用于特定于该实例的操作。我们通常不需要类方法，但如果我们有共享于所有实例的数据，创建类方法比实例方法更有意义。

由于_index_tickers是私有的，我们希望为类的用户提供查看可用项的简便方法。因此，我们将为_index_tickers的键创建一个属性。为此，我们使用@property装饰器。@property和@classmethod）装饰器，并编写我们自己的装饰器以清理并标准化跨方法收集数据的结果（@label_sanitizer）。要使用装饰器，我们将其放在函数或方法定义之上：

    @property
    def available_tickers(self):
        """Indices whose tickers are supported."""
        return list(self._index_tickers.keys())

此外，我们通过类方法提供获取股票代码的方式，因为我们的股票代码存储在类变量中。按照惯例，类方法将cls作为第一个参数，而实例方法将self作为第一个参数：

    @classmethod
    def get_index_ticker(cls, index):
        """
        Get the ticker of the specified index, if known.
        Parameters:
            - index: The name of the index; check 
              `available_tickers` for full list which includes:
                - 'S&P 500' for S&P 500,
                - 'Dow Jones' for Dow Jones Industrial Average,
                - 'NASDAQ' for NASDAQ Composite Index
        Returns: 
            The ticker as a string if known, otherwise `None`.
        """
        try:
            index = index.upper()
        except AttributeError:
            raise ValueError('`index` must be a string')
        return cls._index_tickers.get(index, None)

提示

如果我们想要禁止代码中的某些操作，我们可以检查并根据需要raise错误；这允许我们提供更具信息性的错误消息，或者在重新引发错误之前，简单地附加一些额外的操作（通过raise没有表达式）。如果我们希望在某些事情出错时运行特定代码，则使用try...except块：我们将可能出问题的代码放在try中，并将遇到问题时该做的事情放在except子句中。

当我们进入金融工具的技术分析部分时，我们将需要无风险回报率来计算一些指标。这是指没有金融损失风险的投资回报率；在实际操作中，我们使用 10 年期美国国债收益率。由于这个利率会依赖于我们分析的日期范围，因此我们将把这个功能添加到StockReader类中，避免自己查找。我们将使用pandas_datareader包从圣路易斯联邦储备银行收集数据（fred.stlouisfed.org/series/DGS10），提供选择返回我们研究的日期范围的每日利率（用于分析数据本身），或者仅返回最后一个利率（如果我们需要单一值进行计算）：

    def get_risk_free_rate_of_return(self, last=True):
        """
        Get risk-free rate of return w/ 10-year US T-bill 
        from FRED (https://fred.stlouisfed.org/series/DGS10)
        Parameter:
            - last: If `True`, return the rate on the last
              date in the date range else, return a `Series` 
              object for the rate each day in the date range.
        Returns:
            A single value or a `pandas.Series` object.
        """
        data = web.DataReader(
            'DGS10', 'fred', start=self.start, end=self.end
        )
        data.index.rename('date', inplace=True)
        data = data.squeeze()
        return data.asof(self.end) \
            if last and isinstance(data, pd.Series) else data

剩余的方法代码用pass替代，表示告诉 Python 什么都不做（并提醒我们稍后更新），以便代码可以按原样运行。我们将在下一节编写以下方法：

    @label_sanitizer
    def get_ticker_data(self, ticker):
        pass    
    def get_index_data(self, index):
        pass
    def get_bitcoin_data(self, currency_code):
        pass
    @label_sanitizer 
    def get_forex_rates(self, from_currency, to_currency,
                        **kwargs):
        pass

重要提示

由于我们不会查看外汇汇率，本章不涉及get_forex_rates()方法；然而，这个方法提供了如何使用pandas_datareader包的另一个示例，因此我鼓励你看看它。请注意，为了使用这个方法，你需要从 AlphaVantage 获取一个免费的 API 密钥，网址是www.alphavantage.co/support/#api-key。

get_ticker_data()和get_forex_rates()方法都使用了@label_sanitizer装饰器，这样可以将我们从不同来源收到的数据统一为相同的列名，从而避免我们后续清理数据。@label_sanitizer装饰器定义在stock_analysis/utils.py模块中。和之前一样，我们首先来看一下utils模块的文档字符串和导入：

"""Utility functions for stock analysis."""
from functools import wraps
import re
import pandas as pd

接下来，我们有_sanitize_label()函数，它将清理单个标签。请注意，我们在函数名前加了下划线，因为我们不打算让包的用户直接使用它——它是为我们的装饰器使用的：

def _sanitize_label(label):
    """
    Clean up a label by removing non-letter, non-space 
    characters and putting in all lowercase with underscores
    replacing spaces.
    Parameters:
        - label: The text you want to fix.
    Returns: 
        The sanitized label.
    """
    return re.sub(r'[^\w\s]', '', label)\
        .lower().replace(' ', '_')

最后，我们定义了@label_sanitizer装饰器，它是一个清理我们从互联网上获得的数据的列名和索引名的函数。没有这个装饰器，我们收集的数据中的列名可能会有一些意外的字符，比如星号或空格，使得数据难以使用。通过使用这个装饰器，方法将始终返回一个清理过列名的数据框，省去了我们的一步：

def label_sanitizer(method):
    """
    Decorator around a method that returns a dataframe to
    clean up all labels in said dataframe (column names and 
    index name) by using `_sanitize_label()`.
    Parameters:
        - method: The method to wrap.
    Returns: 
        A decorated method or function.
    """
    @wraps(method) # keep original docstring for help()
    def method_wrapper(self, *args, **kwargs):
        df = method(self, *args, **kwargs)
        # fix the column names
        df.columns = [
            _sanitize_label(col) for col in df.columns
        ]
        # fix the index name
        df.index.rename(
            _sanitize_label(df.index.name), inplace=True
        )
        return df
    return method_wrapper

请注意，label_sanitizer()函数的定义中也有一个装饰器。来自标准库functools模块的@wraps装饰器会将装饰过的函数/方法的文档字符串与原始文档字符串相同；这是必要的，因为装饰操作实际上会创建一个新函数/方法，从而使得help()函数变得无效，除非我们进行干预。

小贴士

使用@label_sanitizer语法的方式是method = label_sanitizer(method)。不过，两者都是有效的。

现在我们已经理解了装饰器，准备好完成StockReader类的构建。请注意，我们还将为stock_analysis包中的其他类使用并创建附加的装饰器，因此在继续之前，请确保你对这些装饰器已经足够熟悉。

从 Yahoo! Finance 收集历史数据

我们的数据收集的基础将是get_ticker_data()方法。它使用pandas_datareader包从 Yahoo! Finance 获取数据：

@label_sanitizer
def get_ticker_data(self, ticker):
    """
    Get historical OHLC data for given date range and ticker.
    Parameter:
        - ticker: The stock symbol to lookup as a string.
    Returns: A `pandas.DataFrame` object with the stock data.
    """
    return web.get_data_yahoo(ticker, self.start, self.end)

重要提示

过去，pandas_datareader和 Yahoo! Finance API 曾出现过一些问题，导致pandas_datareader开发者通过web.DataReader()函数（pandas-datareader.readthedocs.io/en/latest/whatsnew.html#v0-6-0-january-24-2018）弃用了对它的支持；因此，我们必须使用他们的替代方法：web.get_data_yahoo()。

要收集股票市场指数的数据，我们可以使用get_index_data()方法，该方法首先查找指数的股票代码，然后调用我们刚刚定义的get_ticker_data()方法。请注意，由于get_ticker_data()方法使用了@label_sanitizer装饰器，因此get_index_data()方法不需要使用@label_sanitizer装饰器：

def get_index_data(self, index):
    """
    Get historical OHLC data from Yahoo! Finance
    for the chosen index for given date range.
    Parameter:
        - index: String representing the index you want
          data for, supported indices include:
            - 'S&P 500' for S&P 500,
            - 'Dow Jones' for Dow Jones Industrial Average,
            - 'NASDAQ' for NASDAQ Composite Index
    Returns: 
        A `pandas.DataFrame` object with the index data.
    """
    if index not in self.available_tickers:
        raise ValueError(
            'Index not supported. Available tickers'
            f"are: {', '.join(self.available_tickers)}"
        )
    return self.get_ticker_data(self.get_index_ticker(index))

Yahoo! Finance 也提供比特币的数据；然而，我们必须选择一个货币来使用。get_bitcoin_data()方法接受一个货币代码来创建 Yahoo! Finance 搜索的符号（例如，BTC-USD 表示以美元计价的比特币数据）。实际的数据收集仍然由get_ticker_data()方法处理：

def get_bitcoin_data(self, currency_code):
    """
    Get bitcoin historical OHLC data for given date range. 
    Parameter:
        - currency_code: The currency to collect the bitcoin
          data in, e.g. USD or GBP.
    Returns: 
        A `pandas.DataFrame` object with the bitcoin data.
    """
    return self\
        .get_ticker_data(f'BTC-{currency_code}')\
        .loc[self.start:self.end] # clip dates

此时，StockReader类已经可以使用，因此让我们在financial_analysis.ipynb笔记本中开始工作，并导入将用于本章其余部分的stock_analysis包：

>>> import stock_analysis

当我们导入stock_analysis包时，Python 会运行stock_analysis/__init__.py文件：

"""Classes for making technical stock analysis easier."""
from .stock_analyzer import StockAnalyzer, AssetGroupAnalyzer
from .stock_modeler import StockModeler
from .stock_reader import StockReader
from .stock_visualizer import \
    StockVisualizer, AssetGroupVisualizer

重要提示

stock_analysis/__init__.py文件中的代码使我们更容易访问包中的类——例如，我们不需要运行stock_analysis.stock_reader.StockReader()，只需运行stock_analysis.StockReader()即可创建一个StockReader对象。

接下来，我们将通过提供数据收集的开始日期和（可选的）结束日期，创建StockReader类的实例。我们将使用 2019-2020 年的数据。请注意，当我们运行此代码时，Python 会调用StockReader.__init__()方法：

>>> reader = \
...     stock_analysis.StockReader('2019-01-01', '2020-12-31')

现在，我们将收集Facebook、Apple、Amazon、Netflix 和 Google（FAANG）、S&P 500 和比特币数据。由于我们使用的所有股票的价格都是以美元计价的，因此我们会请求以美元计价的比特币数据。请注意，我们使用了生成器表达式和多重赋值来获取每个 FAANG 股票的数据框：

>>> fb, aapl, amzn, nflx, goog = (
...     reader.get_ticker_data(ticker)
...     for ticker in ['FB', 'AAPL', 'AMZN', 'NFLX', 'GOOG']
... )
>>> sp = reader.get_index_data('S&P 500')
>>> bitcoin = reader.get_bitcoin_data('USD') 

提示

确保运行help(stock_analysis.StockReader)或help(reader)，以查看所有已定义的方法和属性。输出会清楚地标明哪些方法是类方法，并将属性列在底部的data descriptors部分。这是熟悉新代码的重要步骤。

探索性数据分析

现在我们有了数据，我们想要熟悉它。正如我们在第五章《使用 Pandas 和 Matplotlib 进行数据可视化》和第六章《使用 Seaborn 绘图与自定义技巧》中看到的，创建好的可视化需要了解matplotlib，并且——根据数据格式和可视化的最终目标——还需要了解seaborn。就像我们在StockReader类中做的那样，我们希望让用户更容易地可视化单个资产和资产组，因此，我们不会指望我们包的用户（以及可能的合作者）精通matplotlib和seaborn，而是将围绕这些功能创建包装器。这意味着用户只需能够使用stock_analysis包来可视化他们的财务数据。此外，我们能够为可视化的外观设定标准，避免在进行每次新分析时复制和粘贴大量代码，从而带来一致性和效率的提升。

为了使这一切成为可能，我们在stock_analysis/stock_visualizer.py中有Visualizer类。这个文件中有三个类：

• Visualizer：这是定义Visualizer对象功能的基类。大多数方法是抽象的，这意味着从这个父类继承的子类（子类）需要重写这些方法并实现代码；这些方法定义了对象应该做什么，但不涉及具体细节。

• StockVisualizer：这是我们将用来可视化单个资产的子类。

• AssetGroupVisualizer：这是我们将用来通过groupby()操作可视化多个资产的子类。

在讨论这些类的代码之前，让我们先来看一下stock_analysis/utils.py文件中的一些附加函数，这些函数将帮助我们创建这些资产组并为 EDA 目的描述它们。对于这些函数，我们需要导入pandas：

import pandas as pd

group_stocks()函数接受一个字典，字典将资产的名称映射到该资产的数据框，并输出一个新的数据框，包含输入数据框的所有数据以及一列新数据，标明数据属于哪个资产：

def group_stocks(mapping):
    """
    Create a new dataframe with many assets and a new column 
    indicating the asset that row's data belongs to.
    Parameters:
        - mapping: A key-value mapping of the form 
                   {asset_name: asset_df}
    Returns: 
        A new `pandas.DataFrame` object
    """
    group_df = pd.DataFrame()
    for stock, stock_data in mapping.items():
        df = stock_data.copy(deep=True)
        df['name'] = stock
        group_df = group_df.append(df, sort=True)
    group_df.index = pd.to_datetime(group_df.index)
    return group_df

由于在整个包中，我们将有许多方法和函数期望它们的数据框具有特定格式，因此我们将构建一个新的装饰器：@validate_df。这个装饰器检查传递给给定方法或函数的输入是否为DataFrame类型的对象，并且至少包含装饰器columns参数指定的列。我们将以set对象的形式提供这些列。这样我们就可以检查我们需要的列与输入数据中的列之间的集合差异（参见第四章，聚合 Pandas 数据框，了解集合操作）。如果数据框包含我们要求的列（至少），那么集合差异将为空，这意味着数据框通过了测试。如果违反了这些条件，装饰器将抛出ValueError。

让我们来看一下在stock_analysis/utils.py文件中是如何定义的：

def validate_df(columns, instance_method=True):
    """
    Decorator that raises a `ValueError` if input isn't a
    `DataFrame` or doesn't contain the proper columns. Note 
    the `DataFrame` must be the first positional argument
    passed to this method.
    Parameters:
        - columns: A set of required column names.
          For example, {'open', 'high', 'low', 'close'}.
        - instance_method: Whether or not the item being
          decorated is an instance method. Pass `False` to 
          decorate static methods and functions.
    Returns:
        A decorated method or function.
    """
    def method_wrapper(method):
        @wraps(method)
        def validate_wrapper(self, *args, **kwargs):
            # functions and static methods don't pass self so
            # self is the 1st positional argument in that case
            df = (self, *args)[0 if not instance_method else 1]
            if not isinstance(df, pd.DataFrame):
                raise ValueError(
                    'Must pass in a pandas `DataFrame`'
                )
            if columns.difference(df.columns):
                raise ValueError(
                    'Dataframe must contain the following'
                    f' columns: {columns}'
                )
            return method(self, *args, **kwargs)
        return validate_wrapper
    return method_wrapper

使用group_stocks()函数创建的组可以通过describe_group()函数在一个输出中进行描述。group_stocks()函数添加了一个名为name的列，而describe_group()会查找这个列，因此我们使用@validate_df装饰器确保格式正确，然后再尝试运行该函数：

@validate_df(columns={'name'}, instance_method=False)
def describe_group(data):
    """
    Run `describe()` on the asset group.
    Parameters:
        - data: Grouped data resulting from `group_stocks()`
    Returns: 
        The transpose of the grouped description statistics.
    """
    return data.groupby('name').describe().T

让我们使用group_stocks()函数为我们的分析创建一些资产组：

>>> from stock_analysis.utils import \
...     group_stocks, describe_group
>>> faang = group_stocks({
...     'Facebook': fb, 'Apple': aapl, 'Amazon': amzn, 
...     'Netflix': nflx, 'Google': goog
... })
>>> faang_sp = group_stocks({
...     'Facebook': fb, 'Apple': aapl, 'Amazon': amzn, 
...     'Netflix': nflx, 'Google': goog, 'S&P 500': sp
... })
>>> all_assets = group_stocks({
...     'Bitcoin': bitcoin, 'S&P 500': sp, 'Facebook': fb, 
...     'Apple': aapl, 'Amazon': amzn, 'Netflix': nflx, 
...     'Google': goog
... })

使用这些组，describe()的输出在比较时要比分别对每个数据框运行它更有信息量。describe_group()函数通过groupby()来运行describe()。这使得查看不同资产的收盘价汇总更加方便：

>>> describe_group(all_assets).loc['close',]

一眼看去，我们可以看到比特币的数据比其他资产更多。这是因为比特币的价格每天都会变化，而股票数据只包含交易日的数据。我们还可以从中得出一个结论，即比特币不仅波动性大，而且价值远高于其他资产：

![图 7.6 – 各金融工具的收盘价汇总统计]

](github.com/OpenDocCN/f…)

图 7.6 – 各金融工具的收盘价汇总统计

如果我们不想单独查看每个资产，我们可以将它们组合成一个投资组合，视其为一个单一资产。stock_analysis/utils.py中的make_portfolio()函数按日期对数据进行分组并对所有列求和，从而得出我们投资组合的总股价和交易量：

@validate_df(columns=set(), instance_method=False)
def make_portfolio(data, date_level='date'):
    """
    Make a portfolio of assets by grouping by date and 
    summing all columns.
    Note: the caller is responsible for making sure the 
    dates line up across assets and handling when they don't.
    """
    return data.groupby(level=date_level).sum()

该函数假设资产是以相同的频率进行交易的。比特币每天交易，而股市则不是。因此，如果我们的投资组合是比特币和股市的混合体，在使用此函数之前，我们需要决定如何处理这一差异；请参考我们在第三章中关于重新索引的讨论，使用 Pandas 进行数据清洗，以寻找可能的策略。我们将在本章末尾的练习中使用此函数，构建一个由 FAANG 股票组成的投资组合，这些股票的交易频率相同，以便观察盘后交易对整个 FAANG 股票的影响。

可视化器类族

正如我们从前几章中了解到的那样，可视化将使我们的分析变得更加轻松，因此让我们开始讨论 stock_analysis/stock_visualizer.py 中的 Visualizer 类。首先，我们将定义我们的基础类 Visualizer。以下的 UML 图告诉我们这是我们的基础类，因为它有箭头指向它，这些箭头来自子类（AssetGroupVisualizer 和 StockVisualizer）：

图 7.7 – 可视化器类层次结构

图 7.7 还告诉我们将为本节中的每个类定义哪些方法。这包括可视化盘后交易影响（after_hours_trades()）和资产价格随时间变化（evolution_over_time()）的方法，我们将用它们来进行资产的可视化比较。

我们以文档字符串和导入开始模块。对于我们的可视化，我们将需要 matplotlib、numpy、pandas 和 seaborn，以及 mplfinance（一个用于金融可视化的 matplotlib 派生包）：

"""Visualize financial instruments."""
import math
import matplotlib.pyplot as plt
import mplfinance as mpf
import numpy as np
import pandas as pd
import seaborn as sns
from .utils import validate_df

接下来，我们开始定义 Visualizer 类。这个类将保存将用于可视化的数据，因此我们将其放入 __init__() 方法中：

class Visualizer:
    """Base visualizer class not intended for direct use."""
    @validate_df(columns={'open', 'high', 'low', 'close'})
    def __init__(self, df):
        """Store the input data as an attribute."""
        self.data = df

该基础类将为我们提供调用 matplotlib 函数所需的功能；静态方法不依赖于类的数据。我们使用 @staticmethod 装饰器定义 add_reference_line() 方法，用于添加水平或垂直线（以及介于两者之间的任何内容）；注意，我们没有将 self 或 cls 作为第一个参数：

    @staticmethod
    def add_reference_line(ax, x=None, y=None, **kwargs):
        """
        Static method for adding reference lines to plots.
        Parameters:
            - ax: `Axes` object to add the reference line to.
            - x, y: The x, y value to draw the line at as a 
              single value or numpy array-like structure.
                - For horizontal: pass only `y`
                - For vertical: pass only `x`
                - For AB line: pass both `x` and `y`
            - kwargs: Additional keyword args. to pass down.
        Returns:
            The matplotlib `Axes` object passed in.
        """
        try:
            # numpy array-like structures -> AB line
            if x.shape and y.shape:
                ax.plot(x, y, **kwargs)
        except:
            # error triggers if x or y isn't array-like
            try:
                if not x and not y:
                    raise ValueError(
                        'You must provide an `x` or a `y`'
                    )
                elif x and not y:
                    ax.axvline(x, **kwargs) # vertical line
                elif not x and y:
                    ax.axhline(y, **kwargs) # horizontal line
            except:
                raise ValueError(
                    'If providing only `x` or `y`, '
                    'it must be a single value'
                )
        ax.legend()
        return ax

提示

有关类方法、静态方法和抽象方法的更多信息，请参见进一步阅读部分。

shade_region() 静态方法用于向图表添加阴影区域，它类似于 add_reference_line() 静态方法：

    @staticmethod
    def shade_region(ax, x=tuple(), y=tuple(), **kwargs):
        """
        Static method for shading a region on a plot.
        Parameters:
            - ax: `Axes` object to add the shaded region to.
            - x: Tuple with the `xmin` and `xmax` bounds for 
              the rectangle drawn vertically.
            - y: Tuple with the `ymin` and `ymax` bounds for 
              the rectangle drawn horizontally.
            - kwargs: Additional keyword args. to pass down.
        Returns: 
            The matplotlib `Axes` object passed in.
        """
        if not x and not y:
            raise ValueError(
                'You must provide an x or a y min/max tuple'
            )
        elif x and y:
            raise ValueError('You can only provide x or y.')
        elif x and not y:
            ax.axvspan(*x, **kwargs) # vertical region
        elif not x and y:
            ax.axhspan(*y, **kwargs) # horizontal region
        return ax

由于我们希望我们的绘图功能具有灵活性，我们将定义一个静态方法，使得我们可以轻松地绘制一个或多个项目，而无需事先检查项目的数量。这将在我们使用 Visualizer 类作为基础构建的类中得到应用：

    @staticmethod
    def _iter_handler(items):
        """
        Static method for making a list out of an item if 
        it isn't a list or tuple already.
        Parameters:
            - items: The variable to make sure it is a list.
        Returns: The input as a list or tuple.
        """
        if not isinstance(items, (list, tuple)):
            items = [items]
        return items

我们希望支持单一资产和资产组的窗口函数；然而，这一实现会有所不同，因此我们将在超类中定义一个抽象方法（一个没有实现的方法），子类将重写它以提供具体实现：

    def _window_calc(self, column, periods, name, func, 
                     named_arg, **kwargs):
        """
        To be implemented by subclasses. Defines how to add 
        lines resulting from window calculations.
        """
        raise NotImplementedError('To be implemented by '
                                  'subclasses.')

这样我们就可以定义依赖于_window_calc()的功能，但不需要知道具体实现，只需要知道结果。moving_average()方法使用_window_calc()将移动平均线添加到图表中：

    def moving_average(self, column, periods, **kwargs):
        """
        Add line(s) for the moving average of a column.
        Parameters:
            - column: The name of the column to plot.
            - periods: The rule or list of rules for 
              resampling, like '20D' for 20-day periods.
            - kwargs: Additional arguments to pass down.
        Returns: A matplotlib `Axes` object.
        """
        return self._window_calc(
            column, periods, name='MA', named_arg='rule',
            func=pd.DataFrame.resample, **kwargs
        )

类似地，我们定义了exp_smoothing()方法，它将使用_window_calc()将指数平滑的移动平均线添加到图表中：

    def exp_smoothing(self, column, periods, **kwargs):
        """
        Add line(s) for the exponentially smoothed moving 
        average of a column.
        Parameters:
            - column: The name of the column to plot.
            - periods: The span or list of spans for,
              smoothing like 20 for 20-day periods.
            - kwargs: Additional arguments to pass down.
        Returns: 
            A matplotlib `Axes` object.
        """
        return self._window_calc(
            column, periods, name='EWMA',
            func=pd.DataFrame.ewm, named_arg='span', **kwargs
        )

请注意，虽然我们有方法可以将移动平均和指数平滑移动平均添加到列的图表中，但它们都调用了_window_calc()，该方法在此处未定义。这是因为每个子类都会有自己的_window_calc()实现，而它们将继承顶层方法，无需重写moving_average()或exp_smoothing()。

重要提示

请记住，以单个下划线（_）开头的方法是 Python 对该类对象的help()版本。我们将_window_calc()作为私有方法创建，因为Visualizer类的用户只需要调用moving_average()和exp_smoothing()。

最后，我们将为所有子类添加占位符方法。这些是抽象方法，将由每个子类单独定义，因为实现会根据我们是可视化单一资产还是资产组而有所不同。为了简洁，以下是该类中定义的一部分抽象方法：

    def evolution_over_time(self, column, **kwargs):
        """Creates line plots."""
        raise NotImplementedError('To be implemented by '
                                  'subclasses.')
    def after_hours_trades(self):
        """Show the effect of after-hours trading."""
        raise NotImplementedError('To be implemented by '
                                  'subclasses.')
    def pairplot(self, **kwargs):
        """Create pairplots."""
        raise NotImplementedError('To be implemented by '
                                  'subclasses.')

子类还将定义它们特有的方法，并/或根据需要重写Visualizer类的实现。它们没有重写的方法将会继承。通过使用Visualizer来定义所有Visualizers应做的事情，然后提供更具体的版本，例如仅处理单一资产的StockVisualizer类。

可视化股票

我们将通过继承Visualizer类来开始实现StockVisualizer类；我们选择不重写__init__()方法，因为StockVisualizer类只会有一个数据框作为属性。相反，我们将为需要添加（该类特有的）或重写的方法提供实现。

重要提示

为了简洁起见，我们只涵盖一部分功能；然而，我强烈建议你阅读完整的代码库并在笔记本中测试功能。

我们将重写的第一个方法是evolution_over_time()，它将创建一个随着时间变化的列的折线图：

class StockVisualizer(Visualizer):
    """Visualizer for a single stock."""
    def evolution_over_time(self, column, **kwargs):
        """
        Visualize the evolution over time of a column.
        Parameters:
            - column: The name of the column to visualize.
            - kwargs: Additional arguments to pass down.
        Returns:
            A matplotlib `Axes` object.
        """
        return self.data.plot.line(y=column, **kwargs)

接下来，我们将使用 mplfinance 创建一个蜡烛图，这是一种将 OHLC 数据一起可视化的方法。每一行 OHLC 时间序列将被绘制为一根蜡烛。当蜡烛图为黑色时，资产的收盘价低于开盘价（表示亏损）；当蜡烛图为白色时，资产的收盘价高于开盘价，如下图所示：

图 7.8 – 理解蜡烛图

candlestick() 方法还提供了重新采样数据、显示交易量和绘制特定日期范围的选项：

    def candlestick(self, date_range=None, resample=None, 
                    volume=False, **kwargs):
        """
        Create a candlestick plot for the OHLC data.
        Parameters:
            - date_range: String or `slice()` of dates to 
              pass to `loc[]`, if `None` the plot will be 
              for the full range of the data.
            - resample: The offset to use for resampling 
              the data, if desired.
            - volume: Whether to show a bar plot for volume 
              traded under the candlesticks
            - kwargs: Keyword args for `mplfinance.plot()`
        """
        if not date_range:
            date_range = slice(
                self.data.index.min(), self.data.index.max()
            )
        plot_data = self.data.loc[date_range]
        if resample:
            agg_dict = {
                'open': 'first', 'close': 'last',
                'high': 'max', 'low': 'min', 'volume': 'sum'
            }
            plot_data = plot_data.resample(resample).agg({
                col: agg_dict[col] for col in plot_data.columns
                if col in agg_dict
            })
        mpf.plot(
            plot_data, type='candle', volume=volume, **kwargs
        )

现在，我们添加 after_hours_trades() 方法，帮助我们可视化盘后交易对单个资产的影响，亏损部分用红色条形图表示，盈利部分用绿色条形图表示：

    def after_hours_trades(self):
        """
        Visualize the effect of after-hours trading.
        Returns: A matplotlib `Axes` object.
        """
        after_hours = self.data.open - self.data.close.shift()
        monthly_effect = after_hours.resample('1M').sum()
        fig, axes = plt.subplots(1, 2, figsize=(15, 3))
        after_hours.plot(
            ax=axes[0],
            title='After-hours trading\n'
                  '(Open Price - Prior Day\'s Close)'
        ).set_ylabel('price')
        monthly_effect.index = \
            monthly_effect.index.strftime('%Y-%b')
        monthly_effect.plot(
            ax=axes[1], kind='bar', rot=90,
            title='After-hours trading monthly effect',
            color=np.where(monthly_effect >= 0, 'g', 'r')
        ).axhline(0, color='black', linewidth=1)
        axes[1].set_ylabel('price')
        return axes

接下来，我们将添加一个静态方法，让我们可以填充两条曲线之间的区域。fill_between() 方法将使用 plt.fill_between() 根据哪条曲线较高来为区域上色，绿色或红色：

    @staticmethod
    def fill_between(y1, y2, title, label_higher, label_lower, 
                     figsize, legend_x):
        """
        Visualize the difference between assets.
        Parameters:
            - y1, y2: Data to plot, filling y2 - y1.
            - title: The title for the plot.
            - label_higher: Label for when y2 > y1.
            - label_lower: Label for when y2 <= y1.
            - figsize: (width, height) for the plot dimensions.
            - legend_x: Where to place legend below the plot.
        Returns: A matplotlib `Axes` object.
        """
        is_higher = y2 - y1 > 0
        fig = plt.figure(figsize=figsize)
        for exclude_mask, color, label in zip(
            (is_higher, np.invert(is_higher)),
            ('g', 'r'),
            (label_higher, label_lower)
        ):
            plt.fill_between(
                y2.index, y2, y1, figure=fig,
                where=exclude_mask, color=color, label=label
            )
        plt.suptitle(title)
        plt.legend(
            bbox_to_anchor=(legend_x, -0.1),
            framealpha=0, ncol=2
        )
        for spine in ['top', 'right']:
            fig.axes[0].spines[spine].set_visible(False)
        return fig.axes[0]

open_to_close() 方法将帮助我们通过 fill_between() 静态方法可视化每日开盘价与收盘价之间的差异。如果收盘价高于开盘价，我们会将区域涂成绿色；如果相反，则涂成红色：

    def open_to_close(self, figsize=(10, 4)):
        """
        Visualize the daily change in price from open to close.
        Parameters:
            - figsize: (width, height) of plot
        Returns:
            A matplotlib `Axes` object.
        """
        ax = self.fill_between(
            self.data.open, self.data.close, 
            figsize=figsize, legend_x=0.67,
            title='Daily price change (open to close)',
            label_higher='price rose', label_lower='price fell'
        )
        ax.set_ylabel('price')
        return ax

除了可视化单个资产的开盘价与收盘价之间的差异外，我们还希望比较不同资产之间的价格。fill_between_other() 方法将帮助我们可视化我们为某个资产创建的可视化工具与另一个资产之间的差异，再次使用 fill_between()。当可视化工具中的资产价格高于另一个资产时，我们会将差异部分标为绿色，低于另一个资产时则标为红色：

    def fill_between_other(self, other_df, figsize=(10, 4)):
        """
        Visualize difference in closing price between assets.
        Parameters:
            - other_df: The other asset's data.
            - figsize: (width, height) for the plot.
        Returns: 
            A matplotlib `Axes` object.
        """
        ax = self.fill_between(
            other_df.open, self.data.close, figsize=figsize, 
            legend_x=0.7, label_higher='asset is higher', 
            label_lower='asset is lower', 
            title='Differential between asset price '
                  '(this - other)'
        )
        ax.set_ylabel('price')
        return ax

终于到了重写 _window_calc() 方法的时候，它定义了如何根据单个资产的窗口计算来添加参考线。注意，我们如何能够使用 pipe() 方法（在 第四章 中介绍，聚合 Pandas DataFrame）使我们的窗口计算图与不同的函数兼容，以及如何使用 _iter_handler() 方法使我们的循环在不检查是否有多个参考线需要绘制的情况下工作：

    def _window_calc(self, column, periods, name, func, 
                     named_arg, **kwargs):
        """
        Helper method for plotting a series and adding
        reference lines using a window calculation.
        Parameters:
            - column: The name of the column to plot.
            - periods: The rule/span or list of them to pass 
              to the resampling/smoothing function, like '20D'
              for 20-day periods (resampling) or 20 for a 
              20-day span (smoothing)
            - name: The name of the window calculation (to 
              show in the legend).
            - func: The window calculation function.
            - named_arg: The name of the argument `periods` 
              is being passed as.
            - kwargs: Additional arguments to pass down.
        Returns:
            A matplotlib `Axes` object.
        """
        ax = self.data.plot(y=column, **kwargs)
        for period in self._iter_handler(periods):
            self.data[column].pipe(
                func, **{named_arg: period}
            ).mean().plot(
                ax=ax, linestyle='--',
                label=f"""{period if isinstance(
                    period, str
                ) else str(period) + 'D'} {name}"""
            )
        plt.legend()
        return ax

到目前为止，每个可视化都涉及单一资产的数据；然而，有时我们希望能够可视化资产之间的关系，因此我们将围绕 seaborn 的 jointplot() 函数构建一个封装器：

    def jointplot(self, other, column, **kwargs):
        """
        Generate a seaborn jointplot for given column in 
        this asset compared to another asset.
        Parameters:
            - other: The other asset's dataframe.
            - column: Column to use for the comparison.
            - kwargs: Keyword arguments to pass down.
        Returns: A seaborn jointplot
        """
        return sns.jointplot(
            x=self.data[column], y=other[column], **kwargs
        )

观察资产之间关系的另一种方式是相关矩阵。DataFrame 对象有一个 corrwith() 方法，可以计算每列与另一个数据框中相同列（按名称）的相关系数。这并不能填充热图所需的矩阵，正如我们在前几章所看到的；它实际上是对角线。correlation_heatmap() 方法创建一个矩阵供 sns.heatmap() 函数使用，并用相关系数填充对角线；然后，它通过遮罩确保仅显示对角线。此外，在计算相关性时，我们将使用每列的每日百分比变化，以处理规模差异（例如，苹果股票价格与亚马逊股票价格之间的差异）：

    def correlation_heatmap(self, other):
        """
        Plot the correlations between this asset and another
        one with a heatmap.
        Parameters:
            - other: The other dataframe.
        Returns: A seaborn heatmap
        """
        corrs = \
            self.data.pct_change().corrwith(other.pct_change())
        corrs = corrs[~pd.isnull(corrs)]
        size = len(corrs)
        matrix = np.zeros((size, size), float)
        for i, corr in zip(range(size), corrs):
            matrix[i][i] = corr
        # create mask to only show diagonal
        mask = np.ones_like(matrix)
        np.fill_diagonal(mask, 0)
        return sns.heatmap(
            matrix, annot=True, center=0, vmin=-1, vmax=1,
            mask=mask, xticklabels=self.data.columns, 
            yticklabels=self.data.columns
        )

现在我们已经了解了 StockVisualizer 类的一些功能，我们可以开始进行探索性分析。让我们创建一个 StockVisualizer 对象，对 Netflix 股票数据进行一些 EDA：

>>> %matplotlib inline
>>> import matplotlib.pyplot as plt
>>> netflix_viz = stock_analysis.StockVisualizer(nflx)

一旦我们用 Netflix 数据框初始化了 StockVisualizer 对象，我们就可以生成许多不同类型的图表。我们不会一一举例说明这个对象可以做什么（这留给你自己去实验），但让我们看看随时间变化的收盘价与一些移动平均线，以研究趋势：

>>> ax = netflix_viz.moving_average('close', ['30D', '90D'])
>>> netflix_viz.shade_region(
...     ax, x=('2019-10-01', '2020-07-01'), 
...     color='blue', alpha=0.1
... )
>>> ax.set(title='Netflix Closing Price', ylabel='price ($)')

这些移动平均线给我们提供了股价曲线的平滑版本。请注意，在阴影区域内，90 日移动平均线像是股票价格的天花板：

](tos-cn-i-73owjymdk6/331929e7261f41a39d173301b055b305)

图 7.9 – Netflix 股票价格与移动平均线

交易者根据手头任务的不同，尝试使用不同周期的移动平均线，例如预测股价上涨（股价上涨）并在股价下跌前做出计划性退出（股价下跌）。其他用途包括通过找到支撑线和阻力线来自动计算支撑位和阻力位（我们在 第六章 ，使用 Seaborn 进行绘图与自定义技巧 中首次看到），通过找到支撑线支撑数据的部分，或找到作为数据天花板的部分。当股价接近支撑位时，价格通常会足够吸引人，促使人们买入，从而推动股价上涨（从支撑位向阻力位移动）。然而，当股价达到阻力位时，通常会促使人们卖出，导致股价下跌（从阻力位远离，向支撑位靠近）。

图 7.10 展示了支撑位（绿色）和阻力位（红色）如何分别作为股票价格的下限和上限；一旦价格触及这些边界之一，它通常会因为股票买卖双方的行动而反弹到相反方向：

](tos-cn-i-73owjymdk6/d13ec4cb1f7045c398fd544b3e7f03f5)

图 7.10 – 2018 年 Netflix 股票的支撑位和阻力位示例

通常，指数加权移动平均线（EWMA）可以提供更好的趋势，因为我们可以对最近的值给予更多的权重。让我们来看看对数据进行指数平滑处理的效果：

>>> ax = netflix_viz.exp_smoothing('close', [30, 90]) 
>>> netflix_viz.shade_region(
...     ax, x=('2020-04-01', '2020-10-01'),
...     color='blue', alpha=0.1
... )
>>> ax.set(title='Netflix Closing Price', ylabel='price ($)')

90 天的 EWMA 看起来在阴影区域内充当了支撑位：

图 7.11 – Netflix 股票价格与 EWMAs

提示

该笔记本包含了一个用于交互式可视化移动平均线和指数加权移动平均线（EWMA）的单元格。我们可以使用这些类型的可视化来确定计算的最佳窗口。请注意，使用此单元格可能需要一些额外的设置，但相关设置已在单元格上方标明。

在 第五章《使用 Pandas 和 Matplotlib 可视化数据》的练习中，我们编写了生成可视化图表的代码，展示了盘后交易对 Facebook 的影响；StockVisualizer 类也具备这个功能。我们使用 after_hours_trades() 方法来查看 Netflix 的盘后交易表现：

>>> netflix_viz.after_hours_trades()

Netflix 在 2019 年第三季度的盘后交易表现不佳：

图 7.12 – 可视化盘后交易对 Netflix 股票的影响

我们可以使用蜡烛图来研究 OHLC 数据。我们将使用 candlestick() 方法为 Netflix 创建一个蜡烛图，并同时绘制交易量的条形图。我们还将把数据重新采样为 2 周的时间间隔，以便更清晰地展示蜡烛图：

>>> netflix_viz.candlestick(
...     resample='2W', volume=True, xrotation=90, 
...     datetime_format='%Y-%b –'
... )

从 图 7.8 中记得，当蜡烛图的实体为白色时，意味着股票价值上涨。请注意，大多数情况下，交易量的尖峰伴随着股票价值的上涨：

图 7.13 – 带有交易量的蜡烛图

提示

交易者使用蜡烛图来寻找并分析资产表现中的模式，这些模式可以帮助做出交易决策。查看这篇文章，了解蜡烛图以及交易者常见的分析模式：www.investopedia.com/trading/candlestick-charting-what-is-it/。

在继续之前，我们需要重置图表的样式。mplfinance 包提供了许多可用的样式选项，因此我们暂时返回到我们熟悉的样式：

>>> import matplotlib as mpl
>>> mpl.rcdefaults()
>>> %matplotlib inline

在之前的章节中，我们已经单独看过一只股票（Facebook），所以我们接下来将从另一个角度进行比较，把 Netflix 和其他股票做对比。我们使用 jointplot() 方法来看 Netflix 与标普 500 的对比：

>>> netflix_viz.jointplot(sp, 'close')

如果我们看一下图表，它们似乎有较弱的正相关关系。在金融分析中，我们可以计算一个叫做beta的指标，表示某个资产与一个指数（例如 S&P 500）之间的相关性。我们将在本章后面的金融工具的技术分析部分计算 beta：

图 7.14 – 将 Netflix 与 S&P 500 进行比较

我们可以使用correlation_heatmap()方法，将 Netflix 和 Amazon 之间的相关性可视化为热力图，使用每列的日百分比变化：

>>> netflix_viz.correlation_heatmap(amzn)

Netflix 和 Amazon 之间存在弱正相关，但仅在 OHLC 数据中：

图 7.15 – Netflix 与 Amazon 的相关性热力图

最后，我们可以使用fill_between_other()方法查看另一项资产与 Netflix 相比在价格上的增长（或下降）。在这里，我们将 Netflix 与特斯拉进行比较，看看一个股票超过另一个股票的例子：

>>> tsla = reader.get_ticker_data('TSLA')
>>> change_date = (tsla.close > nflx.close).idxmax()
>>> ax = netflix_viz.fill_between_other(tsla)
>>> netflix_viz.add_reference_line(
...     ax, x=change_date, color='k', linestyle=':', alpha=0.5,
...     label=f'TSLA > NFLX {change_date:%Y-%m-%d}'
... )

请注意，随着阴影区域接近参考线，其高度逐渐缩小——这表示 Netflix 股票和特斯拉股票之间的差异随着时间的推移而减小。在 2020 年 11 月 11 日，当特斯拉超过 Netflix 时，阴影区域的颜色发生变化（从绿色变为红色），并开始增高，因为特斯拉拉大了差距：

图 7.16 – Netflix 与特斯拉的股票价格差异

到目前为止，我们讨论了如何可视化单一资产——在这种情况下是 Netflix——接下来我们将继续，看看如何使用AssetGroupVisualizer类在资产组之间进行一些 EDA 分析。

可视化多个资产

就像之前那样，我们将从Visualizer类继承并编写我们的文档字符串。请注意，AssetGroupVisualizer类还会跟踪用于groupby()操作的列，因此我们会重写__init__()方法；由于这个更改是对已有代码的补充，因此我们也会调用父类的__init__()方法：

class AssetGroupVisualizer(Visualizer):
    """Visualizes groups of assets in a single dataframe."""
    # override for group visuals
    def __init__(self, df, group_by='name'):
        """This object keeps track of the group by column."""
        super().__init__(df)
        self.group_by = group_by

接下来，我们定义evolution_over_time()方法，以便在单个图中绘制资产组中所有资产的相同列，进行对比分析。由于我们的数据形状不同，这次我们将使用seaborn：

    def evolution_over_time(self, column, **kwargs):
        """
        Visualize the evolution over time for all assets.
        Parameters:
            - column: The name of the column to visualize.
            - kwargs: Additional arguments to pass down.
        Returns: A matplotlib `Axes` object.
        """
        if 'ax' not in kwargs:
            fig, ax = plt.subplots(1, 1, figsize=(10, 4))
        else:
            ax = kwargs.pop('ax')
        return sns.lineplot(
            x=self.data.index, y=column, hue=self.group_by,
            data=self.data, ax=ax, **kwargs
        )

当使用seaborn或仅绘制单个资产时，我们不需要担心子图的布局；然而，对于一些其他资产组的可视化，我们需要一种方法来自动确定合理的子图布局。为此，我们将添加_get_layout()方法，该方法将为给定数量的子图生成所需的Figure和Axes对象（由资产组中唯一资产的数量决定）：

    def _get_layout(self):
        """
        Helper method for getting an autolayout of subplots.
        Returns: `Figure` and `Axes` objects to plot with.
        """
        subplots_needed = self.data[self.group_by].nunique()
        rows = math.ceil(subplots_needed / 2)
        fig, axes = \
            plt.subplots(rows, 2, figsize=(15, 5 * rows))
        if rows > 1:
            axes = axes.flatten()
        if subplots_needed < len(axes):
            # remove excess axes from autolayout
            for i in range(subplots_needed, len(axes)):
                # can't use comprehension here
                fig.delaxes(axes[i])
        return fig, axes

现在，我们需要定义 _window_calc() 如何与资产组一起工作。我们需要使用 _get_layout() 方法为组中的每个资产构建子图：

    def _window_calc(self, column, periods, name, func,  
                     named_arg, **kwargs):
        """
        Helper method for plotting a series and adding
        reference lines using a window calculation.
        Parameters:
            - column: The name of the column to plot.
            - periods: The rule/span or list of them to pass 
              to the resampling/smoothing function, like '20D' 
              for 20-day periods (resampling) or 20 for a 
              20-day span (smoothing)
            - name: The name of the window calculation (to 
              show in the legend).
            - func: The window calculation function.
            - named_arg: The name of the argument `periods` 
              is being passed as.
            - kwargs: Additional arguments to pass down.
        Returns: 
            A matplotlib `Axes` object.
        """
        fig, axes = self._get_layout()
        for ax, asset_name in zip(
            axes, self.data[self.group_by].unique()
        ):
            subset = self.data.query(
                f'{self.group_by} == "{asset_name}"'
            )
            ax = subset.plot(
                y=column, ax=ax, label=asset_name, **kwargs
            )
            for period in self._iter_handler(periods):
                subset[column].pipe(
                    func, **{named_arg: period}
                ).mean().plot(
                    ax=ax, linestyle='--',
                    label=f"""{period if isinstance(
                        period, str
                    ) else str(period) + 'D'} {name}"""
                )
            ax.legend()
         plt.tight_layout()
         return ax

我们可以重写 after_hours_trades() 来可视化盘后交易对资产组的影响，方法是使用子图并遍历组中的资产：

    def after_hours_trades(self):
        """
        Visualize the effect of after-hours trading.
        Returns: A matplotlib `Axes` object.
        """
        num_categories = self.data[self.group_by].nunique()
        fig, axes = plt.subplots(
            num_categories, 2, figsize=(15, 3 * num_categories)
        )
        for ax, (name, data) in zip(
            axes, self.data.groupby(self.group_by)
        ):
            after_hours = data.open - data.close.shift()
            monthly_effect = after_hours.resample('1M').sum()
            after_hours.plot(
                ax=ax[0], 
                title=f'{name} Open Price - Prior Day\'s Close'
            ).set_ylabel('price')
            monthly_effect.index = \
                monthly_effect.index.strftime('%Y-%b')
            monthly_effect.plot(
                ax=ax[1], kind='bar', rot=90,
                color=np.where(monthly_effect >= 0, 'g', 'r'),
                title=f'{name} after-hours trading '
                      'monthly effect'
            ).axhline(0, color='black', linewidth=1)
            ax[1].set_ylabel('price')
        plt.tight_layout()
        return axes

使用 StockVisualizer 类，我们能够生成两只资产收盘价之间的联合图，但在这里我们可以重写 pairplot()，以便查看资产组中资产之间收盘价的关系：

    def pairplot(self, **kwargs):
        """
        Generate a seaborn pairplot for this asset group.
        Parameters:
            - kwargs: Keyword arguments to pass down.
        Returns: A seaborn pairplot
        """
        return sns.pairplot(
            self.data.pivot_table(
                values='close', index=self.data.index, 
                columns=self.group_by
            ), diag_kind='kde', **kwargs
        )

最后，我们添加 heatmap() 方法，它生成所有资产组中收盘价之间相关性的热图：

    def heatmap(self, pct_change=True, **kwargs):
        """
        Generate a heatmap for correlations between assets.
        Parameters:
            - pct_change: Whether to show the correlations 
              of the daily percent change in price.
            - kwargs: Keyword arguments to pass down.
        Returns: A seaborn heatmap
        """
        pivot = self.data.pivot_table(
            values='close', index=self.data.index, 
            columns=self.group_by
        )
        if pct_change:
            pivot = pivot.pct_change()
        return sns.heatmap(
            pivot.corr(), annot=True, center=0, 
            vmin=-1, vmax=1, **kwargs
        )

我们可以使用 heatmap() 方法查看资产之间的日变化百分比对比。这将处理资产间的规模差异（谷歌和亚马逊的股价远高于 Facebook 和苹果，这意味着几美元的涨幅对 Facebook 和苹果来说影响更大）：

>>> all_assets_viz = \
...     stock_analysis.AssetGroupVisualizer(all_assets)
>>> all_assets_viz.heatmap()

苹果与标准普尔 500 指数、Facebook 与谷歌之间有最强的相关性，而比特币与任何资产都没有相关性：

图 7.17 – 资产价格之间的相关性

为了简洁起见，避免展示所有可视化资产组的方法（这些方法会生成大量图形），我将把这部分留给你在笔记本中查看并尝试。不过，让我们结合这些 Visualizers 来看看所有资产随时间的演变：

>>> faang_sp_viz = \
...     stock_analysis.AssetGroupVisualizer(faang_sp)
>>> bitcoin_viz = stock_analysis.StockVisualizer(bitcoin)
>>> fig, axes = plt.subplots(1, 2, figsize=(15, 5))
>>> faang_sp_viz.evolution_over_time(
...     'close', ax=axes[0], style=faang_sp_viz.group_by
... )
>>> bitcoin_viz.evolution_over_time(
...     'close', ax=axes[1], label='Bitcoin'
... )

请注意，比特币在 2020 年底大幅上涨（查看 y 轴的刻度），而亚马逊在 2020 年也经历了显著增长：

图 7.18 – 资产价格随时间的变化

现在我们对数据有了充分的了解，我们可以开始查看一些指标了。请注意，虽然我们只查看并使用了部分代码，我鼓励你在本章的笔记本中尝试所有 Visualizer 类的方法；练习题也将提供额外的机会来使用它们。

金融工具的技术分析

在资产的技术分析中，计算一些指标（如累计回报和波动率）来比较不同资产之间的差异。与本章前两部分一样，我们将编写一个包含类的模块来帮助我们。我们将需要 StockAnalyzer 类来分析单一资产的技术指标，和 AssetGroupAnalyzer 类来分析资产组的技术指标。这些类位于 stock_analysis/stock_analyzer.py 文件中。

与其他模块一样，我们将从文档字符串和导入开始：

"""Classes for technical analysis of assets."""
import math
from .utils import validate_df

StockAnalyzer 类

对于单个资产的分析，我们将构建 StockAnalyzer 类来计算给定资产的指标。下图 UML 显示了它提供的所有指标：

图 7.19 – StockAnalyzer 类的结构

StockAnalyzer 实例将使用我们希望进行技术分析的资产数据进行初始化。这意味着我们的 __init__() 方法需要接受数据作为参数：

class StockAnalyzer:
    """Provides metrics for technical analysis of a stock."""
    @validate_df(columns={'open', 'high', 'low', 'close'})
    def __init__(self, df):
        """Create a `StockAnalyzer` object with OHLC data"""
        self.data = df

我们的大部分技术分析计算将依赖于股票的收盘价，因此，为了避免在所有方法中都写 self.data.close，我们将创建一个属性，使我们能够通过 self.close 访问它。这使得我们的代码更加简洁和易于理解：

    @property
    def close(self):
        """Get the close column of the data."""
        return self.data.close

一些计算还需要 close 列的百分比变化，因此我们将为其创建一个属性，方便访问：

    @property
    def pct_change(self):
        """Get the percent change of the close column."""
        return self.close.pct_change()

由于我们将使用枢轴点来计算支撑和阻力水平，枢轴点是数据中最后一天的最高价、最低价和收盘价的平均值，因此我们也会为它创建一个属性：

    @property
    def pivot_point(self):
        """Calculate the pivot point."""
        return (self.last_close + self.last_high
                + self.last_low) / 3

请注意，我们还使用了其他属性——self.last_close、self.last_high 和 self.last_low——这些属性通过在数据上使用 last() 方法定义，然后选择相应的列并使用 iat[] 获取对应的价格：

    @property
    def last_close(self):
        """Get the value of the last close in the data."""
        return self.data.last('1D').close.iat[0]
    @property
    def last_high(self):
        """Get the value of the last high in the data."""
        return self.data.last('1D').high.iat[0]
    @property
    def last_low(self):
        """Get the value of the last low in the data."""
        return self.data.last('1D').low.iat[0]

现在，我们拥有了计算支撑和阻力所需的一切。我们将在三个不同的水平上计算每个值，其中第一个水平离收盘价最近，第三个水平最远。因此，第一个水平是最具限制性的水平，第三个水平则是最不具限制性的。我们将 resistance() 方法定义如下，允许调用者指定计算的级别：

    def resistance(self, level=1):
        """Calculate the resistance at the given level."""
        if level == 1:
            res = (2 * self.pivot_point) - self.last_low
        elif level == 2:
            res = self.pivot_point \
                  + (self.last_high - self.last_low)
        elif level == 3:
            res = self.last_high \
                  + 2 * (self.pivot_point - self.last_low)
        else:
            raise ValueError('Not a valid level.')
        return res

support() 方法的定义方式类似：

    def support(self, level=1):
        """Calculate the support at the given level."""
        if level == 1:
            sup = (2 * self.pivot_point) - self.last_high
        elif level == 2:
            sup = self.pivot_point \
                  - (self.last_high - self.last_low)
        elif level == 3:
            sup = self.last_low \
                  - 2 * (self.last_high - self.pivot_point)
        else:
            raise ValueError('Not a valid level.')
        return sup

接下来，我们将创建用于分析资产波动性的方法。首先，我们将计算收盘价百分比变化的日标准差，计算时需要指定交易期数。为了确保我们不会使用超出数据中期数的交易期数，我们将定义一个属性，其中包含可以用于此参数的最大值：

    @property
    def _max_periods(self):
        """Get the number of trading periods in the data."""
        return self.data.shape[0]

现在我们已经得到了最大值，我们可以定义 daily_std() 方法，该方法计算每日百分比变化的日标准差：

    def daily_std(self, periods=252):
        """
        Calculate daily standard deviation of percent change.
        Parameters:
            - periods: The number of periods to use for the
              calculation; default is 252 for the trading days 
              in a year. Note if you provide a number greater  
              than the number of trading periods in the data,
              `self._max_periods` will be used instead.
        Returns: The standard deviation
        """
        return self.pct_change\
            [min(periods, self._max_periods) * -1:].std()

虽然 daily_std() 本身很有用，但我们可以更进一步，通过将日标准差乘以一年中交易期数的平方根来计算年化波动性，我们假设一年有 252 个交易日：

    def annualized_volatility(self):
        """Calculate the annualized volatility."""
        return self.daily_std() * math.sqrt(252)

此外，我们可以使用 rolling() 方法来查看滚动波动性：

    def volatility(self, periods=252):
        """Calculate the rolling volatility.
        Parameters:
            - periods: The number of periods to use for the 
              calculation; default is 252 for the trading  
              days in a year. Note if you provide a number  
              greater than the number of trading periods in the
              data, `self._max_periods` will be used instead.
        Returns: A `pandas.Series` object.
        """
        periods = min(periods, self._max_periods)
        return self.close.rolling(periods).std()\
               / math.sqrt(periods)

我们经常需要比较不同资产，因此我们提供了 corr_with() 方法，使用每日百分比变化来计算它们之间的相关性：

    def corr_with(self, other):
        """Calculate the correlations between dataframes.
        Parameters:
            - other: The other dataframe.
        Returns: A `pandas.Series` object
        """
        return \
            self.data.pct_change().corrwith(other.pct_change())

接下来，我们定义一些用于比较资产分散程度的指标。在第一章《数据分析入门》中，我们讨论了变异系数（cv()方法）和分位数离散系数（qcd()方法），我们可以利用这些指标来实现此目标，接下来我们将添加这两者：

    def cv(self):
        """
        Calculate the coefficient of variation for the asset.
        The lower this is, the better the risk/return tradeoff.
        """
        return self.close.std() / self.close.mean()
    def qcd(self):
        """Calculate the quantile coefficient of dispersion."""
        q1, q3 = self.close.quantile([0.25, 0.75])
        return (q3 - q1) / (q3 + q1) 

此外，我们还希望有一种方法来量化资产相对于指数的波动性，例如标准普尔 500 指数（S&P 500），为此我们计算beta()方法，允许用户指定用作基准的指数：

    def beta(self, index):
        """
        Calculate the beta of the asset.
        Parameters:
            - index: The data for the index to compare to.
        Returns: 
            Beta, a float.
        """
        index_change = index.close.pct_change()
        beta = self.pct_change.cov(index_change)\
               / index_change.var()
        return beta

接下来，我们定义一个方法来计算资产的累积回报率，作为一个系列。这被定义为一加上收盘价百分比变化的累积乘积：

    def cumulative_returns(self):
        """Calculate cumulative returns for plotting."""
        return (1 + self.pct_change).cumprod()

接下来我们需要支持的几个指标要求计算投资组合的回报。为了简化问题，我们假设每股没有分红，因此投资组合的回报是从起始价格到结束价格的百分比变化，覆盖的数据时间段内。我们将其定义为静态方法，因为我们需要为一个指数计算该值，而不仅仅是self.data中存储的数据：

    @staticmethod
    def portfolio_return(df):
        """
        Calculate return assuming no distribution per share.
        Parameters:
            - df: The asset's dataframe.
        Returns: The return, as a float.
        """
        start, end = df.close[0], df.close[-1]
        return (end - start) / start

虽然 beta 可以让我们将资产的波动性与指数进行比较，但alpha使我们能够将资产的回报与指数的回报进行比较。为此，我们还需要无风险回报率，即没有财务损失风险的投资的回报率；在实际操作中，我们通常使用美国国债作为参考。计算 alpha 需要计算指数和资产的投资组合回报以及 beta：

    def alpha(self, index, r_f):
        """
        Calculates the asset's alpha.
        Parameters:
            - index: The index to compare to.
            - r_f: The risk-free rate of return.
        Returns: Alpha, as a float.
        """
        r_f /= 100
        r_m = self.portfolio_return(index)
        beta = self.beta(index)
        r = self.portfolio_return(self.data)
        alpha = r - r_f - beta * (r_m - r_f)
        return alpha

提示

r_f by 100 before storing the result back in r_f. It's shorthand for r_f = r_f / 100. Python also has these operators for other arithmetic functions—for example, +=, -=, *=, and %=.

我们还希望添加一些方法，告诉我们资产是否处于熊市或牛市，即在过去 2 个月内，股票价格分别下跌或上涨了 20％以上：

    def is_bear_market(self):
        """
        Determine if a stock is in a bear market, meaning its
        return in the last 2 months is a decline of 20% or more
        """
        return \
            self.portfolio_return(self.data.last('2M')) <= -.2
    def is_bull_market(self):
        """
        Determine if a stock is in a bull market, meaning its
        return in the last 2 months is an increase of >= 20%.
        """
        return \
            self.portfolio_return(self.data.last('2M')) >= .2

最后，我们将添加一个方法来计算夏普比率，该比率告诉我们在承担投资波动性时，相对于无风险回报率，我们所获得的超额回报：

    def sharpe_ratio(self, r_f):
        """
        Calculates the asset's Sharpe ratio.
        Parameters:
            - r_f: The risk-free rate of return.
        Returns: 
            The Sharpe ratio, as a float.
        """
        return (
            self.cumulative_returns().last('1D').iat[0] - r_f
        ) / self.cumulative_returns().std()

花些时间消化本模块中的代码，因为我们将在之前讨论的基础上继续构建。我们不会使用所有这些指标进行技术分析，但我鼓励你在本章的笔记本中尝试这些方法。

AssetGroupAnalyzer类

本节中我们将使用的所有计算都定义在StockAnalyzer类中；然而，为了避免对每个要比较的资产都进行计算，我们还将创建AssetGroupAnalyzer类（在同一个模块中），该类能够为一组资产提供这些指标。

StockAnalyzer和AssetGroupAnalyzer类将共享它们的大部分功能，这为它们的继承设计提供了强有力的论据；然而，有时——如在这种情况下——组合设计更有意义。当对象包含其他类的实例时，这就是AssetGroupAnalyzer类的情况：

图 7.20 – AssetGroupAnalyzer 类的结构

我们通过提供资产的数据框和分组列的名称（如果不是 name）来创建 AssetGroupAnalyzer 实例。初始化时，会调用 _composition_handler() 方法来创建 StockAnalyzer 对象的字典（每个资产一个）：

class AssetGroupAnalyzer:
    """Analyzes many assets in a dataframe."""
    @validate_df(columns={'open', 'high', 'low', 'close'})
    def __init__(self, df, group_by='name'):
        """
        Create an `AssetGroupAnalyzer` object with a 
        dataframe of OHLC data and column to group by.
        """
        self.data = df 
        if group_by not in self.data.columns:
            raise ValueError(
                f'`group_by` column "{group_by}" not in df.'
) 
        self.group_by = group_by
        self.analyzers = self._composition_handler()
    def _composition_handler(self):
        """
        Create a dictionary mapping each group to its analyzer,
        taking advantage of composition instead of inheritance.
        """
        return {
            group: StockAnalyzer(data)
            for group, data in self.data.groupby(self.group_by)
        }

AssetGroupAnalyzer 类只有一个公共方法，analyze()—所有实际的计算都委托给存储在 analyzers 属性中的 StockAnalyzer 对象：

    def analyze(self, func_name, **kwargs):
        """
        Run a `StockAnalyzer` method on all assets.
        Parameters:
            - func_name: The name of the method to run.
            - kwargs: Additional arguments to pass down.
        Returns: 
            A dictionary mapping each asset to the result 
            of the calculation of that function.
        """
        if not hasattr(StockAnalyzer, func_name):
            raise ValueError(
                f'StockAnalyzer has no "{func_name}" method.'
            )
        if not kwargs:
            kwargs = {}
        return {
            group: getattr(analyzer, func_name)(**kwargs)
            for group, analyzer in self.analyzers.items()
        }

使用继承时，在这种情况下，所有方法都必须被重写，因为它们无法处理 groupby() 操作。相反，使用组合时，只需要为每个资产创建 StockAnalyzer 对象，并使用字典推导式来进行计算。另一个很棒的地方是，使用 getattr() 时，无需在 AssetGroupAnalyzer 类中镜像方法，因为 analyze() 可以通过名称使用 StockAnalyzer 对象来抓取方法。

比较资产

让我们使用 AssetGroupAnalyzer 类来比较我们收集的所有资产数据。与之前的章节一样，我们不会在这里使用 StockAnalyzer 类中的所有方法，因此请务必自己尝试：

>>> all_assets_analyzer = \
...     stock_analysis.AssetGroupAnalyzer(all_assets)

记住在第一章《数据分析导论》中提到，变异系数（CV）是标准差与均值的比率；这有助于我们比较资产收盘价的变化程度，即使它们的均值差异较大（例如，亚马逊和苹果）。CV 还可以用来比较投资的波动性与预期回报，并量化风险与回报的权衡。让我们使用 CV 来查看哪种资产的收盘价波动最大：

>>> all_assets_analyzer.analyze('cv')
{'Amazon': 0.2658012522278963,
 'Apple': 0.36991905161737615,
 'Bitcoin': 0.43597652683008137,
 'Facebook': 0.19056336194852783,
 'Google': 0.15038618497328074,
 'Netflix': 0.20344854330432688,
 'S&P 500': 0.09536374658108937}

比特币有着最广泛的价格波动，这应该不是什么惊讶的事。与其使用收盘价，不如使用每日的百分比变化来计算年化波动率。这涉及到计算过去一年中百分比变化的标准差，并将其乘以一年中交易日数量的平方根（代码假设为 252）。通过使用百分比变化，相对于资产价格的较大价格波动会受到更严厉的惩罚。使用年化波动率时，Facebook 看起来比我们使用 CV 时波动性更大（尽管仍然不是最波动的资产）：

>>> all_assets_analyzer.analyze('annualized_volatility')
{'Amazon': 0.3851099077041784,
 'Apple': 0.4670809643500882,
 'Bitcoin': 0.4635140114227397,
 'Facebook': 0.45943066572169544,
 'Google': 0.3833720603377728,
 'Netflix': 0.4626772090887299,
 'S&P 500': 0.34491195196047003}

鉴于所有资产在数据集末尾都获得了增值，接下来让我们检查这些资产是否进入了牛市，即在过去 2 个月内，资产的回报增长达到了 20% 或更高：

>>> all_assets_analyzer.analyze('is_bull_market')
{'Amazon': False,
 'Apple': True,
 'Bitcoin': True,
 'Facebook': False,
 'Google': False,
 'Netflix': False,
 'S&P 500': False}

看起来苹果和比特币在 2020 年 11 月和 12 月表现相当突出。其他资产的表现似乎不太好；然而，它们都没有进入熊市（我们可以通过将 'is_bear_market' 传递给 analyze() 来确认这一点）。另一种分析波动性的方法是通过计算贝塔值来将资产与指数进行比较。大于 1 的正值表示波动性高于该指数，而小于-1 的负值则表示与该指数的反向关系：

>>> all_assets_analyzer.analyze('beta', index=sp)
{'Amazon': 0.7563691182389207,
 'Apple': 1.173273501105916,
 'Bitcoin': 0.3716024282483362,
 'Facebook': 1.024592821854751,
 'Google': 0.98620762504024,
 'Netflix': 0.7408228073823271,
 'S&P 500': 1.0000000000000002}

使用之前结果中的贝塔值，我们可以看到与标准普尔 500 指数相比，苹果的波动性最大，这意味着如果这是我们的投资组合（暂时不考虑比特币），添加苹果股票会增加投资组合的风险。然而，我们知道比特币与标准普尔 500 指数之间没有相关性（请参见图 7.17中的相关性热图），所以这个低贝塔值是具有误导性的。

我们将要查看的最后一个指标是 r_f）；我们通常使用美国国债的回报率作为这个数字。利率可以在www.treasury.gov/resource-center/data-chart-center/interest-rates/pages/TextView.aspx?data=yield上查找；或者，我们可以使用 StockReader 对象（reader）来为我们收集这个数据。让我们使用标准普尔 500 指数作为基准，比较这些资产的阿尔法值：

>>> r_f = reader.get_risk_free_rate_of_return() # 0.93
>>> all_assets_analyzer.analyze('alpha', index=sp, r_f=r_f)
{'Amazon': 0.7383391908270172,
 'Apple': 1.7801122522388666,
 'Bitcoin': 6.355297988074054,
 'Facebook': 0.5048625273190841,
 'Google': 0.18537197824248092,
 'Netflix': 0.6500392764754642,
 'S&P 500': -1.1102230246251565e-16}

所有资产的表现都超过了标准普尔 500 指数，而标准普尔 500 作为一个由 500 只股票组成的投资组合，其风险较低、回报较低，这是由于Cycler对象的存在（matplotlib.org/cycler/)，该对象改变颜色和线条样式：

>>> from cycler import cycler
>>> bw_viz_cycler = (
...     cycler(color=[plt.get_cmap('tab10')(x/10)
...                   for x in range(10)])
...     + cycler(linestyle=['dashed', 'solid', 'dashdot',
...                         'dotted', 'solid'] * 2))
>>> fig, axes = plt.subplots(1, 2, figsize=(15, 5))
>>> axes[0].set_prop_cycle(bw_viz_cycler)
>>> cumulative_returns = \
...     all_assets_analyzer.analyze('cumulative_returns')
>>> for name, data in cumulative_returns.items():
...     data.plot(
...         ax=axes[1] if name == 'Bitcoin' else axes[0], 
...         label=name, legend=True
...     )
>>> fig.suptitle('Cumulative Returns')

尽管 2020 年初面临困境，但所有资产都获得了价值。请注意，比特币子图的 y 轴范围是从 0 到 7（右侧子图），而股市子图（左侧）的范围则覆盖了该范围的一半：

图 7.21 – 所有资产的累计回报

现在我们已经对如何分析金融工具有了良好的理解，接下来我们来尝试预测未来的表现。

使用历史数据建模表现

本节的目标是让我们了解如何构建一些模型；因此，以下示例并不意味着是最佳模型，而是为了学习目的的简单且相对快速的实现。再次提醒，stock_analysis 包中有一个适用于本节任务的类：StockModeler。

重要提示

要完全理解本节的统计元素和建模一般原理，我们需要对统计学有一个扎实的理解；然而，本讨论的目的是展示如何将建模技术应用于金融数据，而不深入探讨背后的数学原理。

StockModeler 类

StockModeler 类将使我们更容易构建和评估一些简单的金融模型，而无需直接与 statsmodels 包进行交互。此外，我们将减少生成模型所需的步骤。以下 UML 图显示这是一个相当简单的类。请注意，我们没有属性，因为 StockModeler 是一个静态类（意味着我们不实例化它）：

图 7.22 – StockModeler 类的结构

StockModeler 类定义在 stock_analysis/stock_modeler.py 中，并具有构建模型和进行一些初步性能分析的方法。像往常一样，我们以文档字符串和导入开始模块：

"""Simple time series modeling for stocks."""
import matplotlib.pyplot as plt
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.seasonal import seasonal_decompose
import statsmodels.api as sm
from .utils import validate_df

接下来，我们将开始 StockModeler 类，并在有人尝试实例化它时引发错误：

class StockModeler:
    """Static methods for modeling stocks."""
    def __init__(self):
        raise NotImplementedError(
            "This class must be used statically: " 
            "don't instantiate it."
        )

我们希望这个类支持的任务之一是时间序列分解，我们在第一章《数据分析简介》中讨论过这个内容。我们从 statsmodels 导入了 seasonal_decompose() 函数，所以我们只需在 decompose() 方法中对收盘价调用它：

    @staticmethod
    @validate_df(columns={'close'}, instance_method=False)
    def decompose(df, period, model='additive'):
        """
        Decompose the closing price of the stock into 
        trend, seasonal, and remainder components.
        Parameters:
            - df: The dataframe containing the stock closing
              price as `close` and with a time index.
            - period: The number of periods in the frequency.
            - model: How to compute the decomposition
              ('additive' or 'multiplicative')
        Returns:
            A `statsmodels` decomposition object.
        """
        return seasonal_decompose(
            df.close, model=model, period=period
        )

请注意，我们对 decompose() 方法有两个装饰器。最上面的装饰器应用于它下面装饰器的结果。在这个例子中，我们有以下内容：

staticmethod(
    validate_df(
        decompose, columns={'close'}, instance_method=False
    )
)

此外，我们还希望支持创建 ARIMA 模型，我们在第一章《数据分析简介》中也讨论过这一点。ARIMA 模型使用 ARIMA(p, d, q) 符号，其中 p 是 AR 模型的时间滞后数（或阶数），d 是从数据中减去的过去值的数量（即 I 模型），q 是 MA 模型中使用的周期数。因此，ARIMA(1, 1, 1) 是一个包含一个自回归部分时间滞后、数据差分一次以及一个 1 期的移动平均的模型。如果我们有任何阶数为零的情况，可以去掉它们——例如，ARIMA(1, 0, 1) 等同于 ARMA(1, 1)，而 ARIMA(0, 0, 3) 等同于 MA(3)。季节性 ARIMA 模型写作 *ARIMA(p, d, q)(P, D, Q)*m，其中 m 是季节性模型中的周期数，P、D 和 Q 是季节性 ARIMA 模型的阶数。StockModeler.arima() 方法不支持季节性组件（为了简化）并且接收 p、d 和 q 作为参数，但为了避免混淆，我们会根据它们所代表的 ARIMA 特性给它们命名——例如，ar 代表自回归（p）。此外，我们还将使我们的静态方法在返回之前提供拟合模型的选项：

    @staticmethod
    @validate_df(columns={'close'}, instance_method=False)
    def arima(df, *, ar, i, ma, fit=True, freq='B'):
        """
        Create an ARIMA object for modeling time series.
        Parameters:
            - df: The dataframe containing the stock closing
              price as `close` and with a time index.
            - ar: The autoregressive order (p).
            - i: The differenced order (q).
            - ma: The moving average order (d).
            - fit: Whether to return the fitted model
            - freq: Frequency of the time series
        Returns: 
            A `statsmodels` ARIMA object which you can use 
            to fit and predict.
        """
        arima_model = ARIMA(
            df.close.asfreq(freq).fillna(method='ffill'), 
            order=(ar, i, ma)
        )
        return arima_model.fit() if fit else arima_model

提示

注意方法签名（df, *, ar, i, ma, ...）中有一个星号（*）。这强制要求在调用方法时，列出的参数必须作为关键字参数传递。这是一种确保使用者明确表达他们需求的好方式。

为了配合这一点，我们需要一种评估 ARIMA 模型预测结果的方法，因此我们将添加 arima_predictions() 静态方法。我们还将提供将预测结果返回为 Series 对象或图表的选项：

    @staticmethod
    @validate_df(columns={'close'}, instance_method=False)
    def arima_predictions(df, arima_model_fitted, start, end,       
                          plot=True, **kwargs):
        """
        Get ARIMA predictions as a `Series` object or plot.
        Parameters:
            - df: The dataframe for the stock.
            - arima_model_fitted: The fitted ARIMA model.
            - start: The start date for the predictions.
            - end: The end date for the predictions.
            - plot: Whether to plot the result, default is
              `True` meaning the plot is returned instead of
              the `Series` object containing the predictions.
            - kwargs: Additional arguments to pass down.
        Returns: 
            A matplotlib `Axes` object or predictions 
            depending on the value of the `plot` argument.
        """
        predictions = \
            arima_model_fitted.predict(start=start, end=end)
        if plot:
            ax = df.close.plot(**kwargs)
            predictions.plot(
                ax=ax, style='r:', label='arima predictions'
            )
            ax.legend()
        return ax if plot else predictions

类似于我们为 ARIMA 模型构建的内容，我们还将提供 regression() 方法，用于构建以 1 的滞后期为基础的收盘价线性回归模型。为此，我们将再次使用 statsmodels（在 第九章《在 Python 中入门机器学习》中，我们将使用 scikit-learn 进行线性回归）：

    @staticmethod
    @validate_df(columns={'close'}, instance_method=False)
    def regression(df):
        """
        Create linear regression of time series with lag=1.
        Parameters:
            - df: The dataframe with the stock data.
        Returns: 
            X, Y, and the fitted model
        """
        X = df.close.shift().dropna()
        Y = df.close[1:]
        return X, Y, sm.OLS(Y, X).fit()

与 arima_predictions() 方法类似，我们希望提供一种方式来查看模型的预测结果，可以选择以 Series 对象或图表的形式呈现。与 ARIMA 模型不同，它每次只预测一个值。因此，我们将从最后一个收盘价后的第二天开始进行预测，并通过迭代使用前一个预测值来预测下一个值。为了处理这一切，我们将编写 regression_predictions() 方法：

    @staticmethod
    @validate_df(columns={'close'}, instance_method=False)
    def regression_predictions(df, model, start, end, 
                               plot=True, **kwargs):
        """
        Get linear regression predictions as a `pandas.Series`
        object or plot.
        Parameters:
            - df: The dataframe for the stock.
            - model: The fitted linear regression model.
            - start: The start date for the predictions.
            - end: The end date for the predictions.
            - plot: Whether to plot the result, default is
              `True` meaning the plot is returned instead of
              the `Series` object containing the predictions.
            - kwargs: Additional arguments to pass down.
        Returns: 
            A matplotlib `Axes` object or predictions 
            depending on the value of the `plot` argument.
        """
        predictions = pd.Series(
            index=pd.date_range(start, end), name='close'
        )
        last = df.last('1D').close
        for i, date in enumerate(predictions.index):
            if not i:
                pred = model.predict(last)
            else:
                pred = model.predict(predictions.iloc[i - 1])
            predictions.loc[date] = pred[0]
        if plot:
            ax = df.close.plot(**kwargs)
            predictions.plot(
                ax=ax, style='r:', 
                label='regression predictions'
            )
            ax.legend()
        return ax if plot else predictions

最后，对于 ARIMA 和线性回归模型，我们都希望可视化预测中的误差，或者 resid 属性，它将给出残差；我们只需要将它们绘制为散点图来检查其方差，并使用 KDE 检查其均值。为此，我们将添加 plot_residuals() 方法：

    @staticmethod
    def plot_residuals(model_fitted, freq='B'):
        """
        Visualize the residuals from the model.
        Parameters:
            - model_fitted: The fitted model
            - freq: Frequency that the predictions were 
              made on. Default is 'B' (business day).
        Returns: 
            A matplotlib `Axes` object.
        """
        fig, axes = plt.subplots(1, 2, figsize=(15, 5))
        residuals = pd.Series(
            model_fitted.resid.asfreq(freq), name='residuals'
        )
        residuals.plot(
            style='bo', ax=axes[0], title='Residuals'
        )
        axes[0].set(xlabel='Date', ylabel='Residual')
        residuals.plot(
            kind='kde', ax=axes[1], title='Residuals KDE'
        )
        axes[1].set_xlabel('Residual')
        return axes

现在，让我们使用 Netflix 数据再次体验一下 StockModeler 类。

时间序列分解

如在 第一章《数据分析简介》中提到的，时间序列可以利用指定的频率分解为趋势、季节性和剩余成分。这可以通过 statsmodels 包实现，StockModeler.decompose() 正在使用该包：

>>> from stock_analysis import StockModeler
>>> decomposition = StockModeler.decompose(nflx, 20)
>>> fig = decomposition.plot()
>>> fig.suptitle(
...     'Netflix Stock Price Time Series Decomposition', y=1
... )
>>> fig.set_figheight(6)
>>> fig.set_figwidth(10)
>>> fig.tight_layout()

这将返回 Netflix 股票价格的分解图，频率为 20 个交易日：

图 7.23 – Netflix 股票价格随时间变化的时间序列分解

对于更复杂的模型，我们可以先进行分解，然后围绕这些组成部分构建我们的模型。然而，这超出了本章的范围，因此我们将继续讨论 ARIMA 模型。

ARIMA

正如我们在第一章《数据分析导论》中讨论的那样，ARIMA 模型包含自回归、差分和滑动平均部分。它们也可以使用statsmodels包来构建，StockModeler.arima()方法正是利用这个包；此方法会根据提供的规格返回一个拟合的 ARIMA 模型。这里，我们将使用%%capture魔法命令来避免显示 ARIMA 模型拟合过程中产生的任何警告，因为我们正在构建一个简单的模型来探索其功能：

>>> %%capture
>>> arima_model = StockModeler.arima(nflx, ar=10, i=1, ma=5)

提示

我们选择这些值是因为它们在合理的时间内运行。实际上，我们可以使用在第五章《使用 Pandas 和 Matplotlib 进行数据可视化》中介绍的pandas.plotting模块中的autocorrelation_plot()函数，帮助找到一个合适的ar值。

一旦模型拟合完成，我们可以通过模型的summary()方法获取相关信息：

>>> print(arima_model.summary())

该摘要相当全面，我们在解释时应参考文档；然而，本文可能是更易理解的入门介绍：medium.com/analytics-vidhya/interpreting-arma-model-results-in-statsmodels-for-absolute-beginners-a4d22253ad1c。需要注意的是，解释此摘要需要扎实的统计学基础：

图 7.24 – ARIMA 模型的摘要

为了分析这个模型，我们可以用一种更简单的方法，通过查看StockModeler.plot_residuals()方法来帮助我们可视化地检查残差：

>>> StockModeler.plot_residuals(arima_model)

尽管残差的均值为 0（右侧子图），它们是异方差的——注意到它们的方差随时间增加（左侧子图）：

图 7.25 – 评估我们 ARIMA 模型的残差

提示

当我们查看图 7.24中的模型摘要时，statsmodels使用默认显著性水平 0.05 进行了异方差性的统计检验。检验统计量标记为异方差性（H），p 值标记为Prob(H)（双侧）。需要注意的是，结果在统计上是显著的（p 值小于或等于显著性水平），这意味着我们的残差不太可能是同方差的。

作为构建 ARIMA 模型的替代方案，StockModeler类还提供了使用线性回归来模拟金融工具收盘价的选项。

使用 statsmodels 进行线性回归

StockModeler.regression()方法使用statsmodels构建一个以前一天收盘价为自变量的线性回归模型，预测收盘价：

>>> X, Y, lm = StockModeler.regression(nflx)
>>> print(lm.summary())

再次，summary()方法会给我们模型拟合的统计信息：

图 7.26 – 我们的线性回归模型总结

提示

查看这篇文章以获取如何解读总结的指导：medium.com/swlh/interpreting-linear-regression-through-statsmodels-summary-4796d359035a。

调整后的 R²使得这个模型看起来非常好，因为它接近 1（在第九章，《Python 机器学习入门》中，我们将进一步讨论这个指标）；然而，我们知道这只是因为股票数据高度自相关，所以让我们再次查看残差：

>>> StockModeler.plot_residuals(lm)

这个模型也存在异方差性问题：

图 7.27 – 评估我们线性回归模型的残差

现在让我们看看 ARIMA 模型和线性回归模型在预测 Netflix 股票收盘价方面，哪一个表现更好。

比较模型

为了比较我们的模型，我们需要在一些新的数据上测试它们的预测。让我们收集 2021 年 1 月前两周 Netflix 股票的每日收盘价，并使用StockModeler类中的预测方法来可视化我们的模型预测与现实的对比：

>>> import datetime as dt
>>> start = dt.date(2021, 1, 1)
>>> end = dt.date(2021, 1, 14)
>>> jan = stock_analysis.StockReader(start, end)\
...     .get_ticker_data('NFLX')
>>> fig, axes = plt.subplots(1, 2, figsize=(15, 5))
>>> arima_ax = StockModeler.arima_predictions(
...     nflx, arima_model, start=start, end=end, 
...     ax=axes[0], title='ARIMA', color='b'
... )
>>> jan.close.plot(
...     ax=arima_ax, style='b--', label='actual close'
... )
>>> arima_ax.legend()
>>> arima_ax.set_ylabel('price ($)')
>>> linear_reg = StockModeler.regression_predictions(
...     nflx, lm, start=start, end=end,
...     ax=axes[1], title='Linear Regression', color='b'
... )
>>> jan.close.plot(
...     ax=linear_reg, style='b--', label='actual close'
... )
>>> linear_reg.legend()
>>> linear_reg.set_ylabel('price ($)')

ARIMA 模型的预测看起来更符合我们预期的模式，但鉴于股市的不可预测性，两个模型都远离了 2021 年 1 月前两周实际发生的情况：

图 7.28 – 模型预测与现实对比

正如我们所看到的，预测股票表现并不容易，即使是短短几天的预测也很困难。有许多数据没有被这些模型捕捉到，比如新闻报道、法规以及管理层的变化等。无论模型看起来拟合得多好，都要小心相信这些预测，因为它们仅仅是外推，而且有很多随机性没有被考虑在内。

为了进一步说明这一点，请看以下通过随机游走和股票数据生成的一组图表。只有一组是真实数据，但是哪一组呢？答案会在图表之后出现，所以在查看之前一定要猜一下：

图 7.29 – 真实还是伪造的股票数据？

这些时间序列都来源于同一个点（2019 年 7 月 1 日微软的收盘价），但只有A是真实的股票数据——B、C和D都是随机游走。很难（或者不可能）分辨，对吧？

总结

在这一章中，我们看到构建 Python 包来进行分析应用程序的开发，可以让其他人非常容易地进行他们自己的分析并复现我们的结果，同时也让我们为未来的分析创建可重复的工作流程。

我们在本章创建的stock_analysis包包含了多个类，用于从互联网收集股票数据（StockReader）；可视化单个资产或它们的组合（Visualizer家族）；计算单个资产或它们组合的指标以进行比较（分别为StockAnalyzer和AssetGroupAnalyzer）；以及使用分解、ARIMA 和线性回归进行时间序列建模（StockModeler）。我们还首次了解了如何在StockModeler类中使用statsmodels包。本章展示了我们迄今为止在本书中学习的pandas、matplotlib、seaborn和numpy功能如何结合在一起，以及这些库如何与其他包协同工作以实现自定义应用。我强烈建议你重新阅读stock_analysis包中的代码，并测试一些我们在本章中未涵盖的方法，以确保你掌握了相关概念。

在下一章中，我们将进行另一个应用案例，学习如何构建一个登录尝试的模拟器，并尝试基于规则的异常检测。

练习

使用stock_analysis包完成以下练习。除非另有说明，数据应使用 2019 年到 2020 年末的数据。如果使用StockReader类收集数据时遇到问题，可以在exercises/目录中找到备份的 CSV 文件：

1. 使用StockAnalyzer和StockVisualizer类，计算并绘制 Netflix 收盘价的三个支撑位和阻力位。

2. 使用StockVisualizer类，查看盘后交易对 FAANG 股票的影响：

   a) 作为单独的股票

   b) 使用来自stock_analysis.utils模块的make_portfolio()函数作为组合

3. 使用StockVisualizer.open_to_close()方法，创建一个图表，显示 FAANG 股票每个交易日的开盘价和收盘价之间的区域：如果价格下降，用红色填充；如果价格上涨，用绿色填充。作为附加任务，对比特币和 S&P 500 的组合也做同样的操作。

4. 共同基金和AssetGroupAnalyzer类。

5. 编写一个函数，返回一个包含以下列的单行数据框：alpha、beta、sharpe_ratio、annualized_volatility、is_bear_market和is_bull_market，这些列分别包含通过StockAnalyzer类对给定股票运行相应方法的结果。在AssetGroupAnalyzer.analyze()方法中使用的字典推导式和getattr()函数将会很有用。

6. 使用StockModeler类，在 2019 年 1 月 1 日至 2020 年 11 月 30 日的 S&P 500 数据上构建 ARIMA 模型，并使用该模型预测 2020 年 12 月的表现。务必检查残差，并将预测表现与实际表现进行比较。

7. 请求 AlphaVantage 的 API 密钥 (www.alphavantage.co/support/#api-key)，并使用你之前为收集数据创建的 StockReader 对象，通过 get_forex_rates() 方法收集从 USD 到 JPY 的每日外汇汇率。使用 2019 年 2 月到 2020 年 1 月的数据，构建一个蜡烛图，并将数据重采样为每周一次。提示：查看标准库中的 slice() 函数（docs.python.org/3/library/functions.html#slice），以提供日期范围。

深入阅读

查看以下资源，以获取更多关于本章内容的信息：

• Python 函数装饰器指南: www.thecodeship.com/patterns/guide-to-python-function-decorators/

• Alpha: www.investopedia.com/terms/a/alpha.asp

• Python 中的类和继承入门: www.jesshamrick.com/2011/05/18/an-introduction-to-classes-and-inheritance-in-python/

• Beta: www.investopedia.com/terms/b/beta.asp

• 变异系数（CV）: www.investopedia.com/terms/c/coefficientofvariation.asp

• 类（Python 文档）: docs.python.org/3/tutorial/classes.html

• 盘后交易如何影响股票价格: www.investopedia.com/ask/answers/05/saleafterhours.asp

• 如何创建一个 Python 包: www.pythoncentral.io/how-to-create-a-python-package/

• 如何在 Python 中创建 ARIMA 模型进行时间序列预测: machinelearningmastery.com/arima-for-time-series-forecasting-with-python/

• 使用 statsmodels 的线性回归（Python）: datatofish.com/statsmodels-linear-regression/

• 面向对象编程: python.swaroopch.com/oop.html

• 随机游走: en.wikipedia.org/wiki/Random_walk

• 股票分析: www.investopedia.com/terms/s/stock-analysis.asp

• 支撑和阻力基础知识: www.investopedia.com/trading/support-and-resistance-basics/

• 技术分析: www.investopedia.com/technical-analysis-4689657

• 如何在 Python 中使用静态方法、类方法或抽象方法的权威指南: julien.danjou.info/guide-python-static-class-abstract-methods/

• 编写设置脚本: docs.python.org/3/distutils/setupscript.html

第九章：第八章：基于规则的异常检测

是时候抓住一些试图通过暴力破解攻击访问网站的黑客了——他们通过尝试一堆用户名和密码组合，直到成功登录为止。这种攻击非常嘈杂，因此为我们提供了大量的数据点用于异常检测，即寻找由与我们认为典型活动不同的过程产生的数据。黑客将被模拟，并不会像现实中那样狡猾，但这将为我们提供良好的异常检测实践机会。

我们将创建一个包，处理登录尝试的模拟，以便生成本章的数据。知道如何进行模拟是我们工具箱中一个至关重要的技能。有时，用精确的数学解法解决问题很困难；然而，定义系统中小组件如何工作可能比较容易。在这种情况下，我们可以模拟小组件并模拟整个系统的行为。模拟结果为我们提供了一个近似解，这可能足以满足我们的目的。

我们将利用基于规则的异常检测来识别模拟数据中的可疑活动。在本章结束时，我们将了解如何使用从各种概率分布生成的随机数来模拟数据，进一步了解 Python 标准库，积累更多构建 Python 包的经验，练习进行探索性数据分析，并初步接触异常检测。

本章将涵盖以下主题：

• 模拟登录尝试以创建本章的数据集

• 执行探索性数据分析以理解模拟数据

• 使用规则和基准线进行异常检测

本章材料

我们将构建一个模拟包来生成本章的数据；该包在 GitHub 上的地址是 github.com/stefmolin/login-attempt-simulator/tree/2nd_edition。在我们设置环境时，这个包已从 GitHub 安装，相关内容在第一章，数据分析导论中有提及；然而，你也可以按照第七章，金融分析——比特币与股市中的说明，安装一个可以编辑的版本。

本章节的代码库可以在github.com/stefmolin/Hands-On-Data-Analysis-with-Pandas-2nd-edition/tree/master/ch_08找到，包含了我们实际分析中使用的笔记本（anomaly_detection.ipynb）、我们将在logs/文件夹中处理的数据文件、用于模拟的数据（位于user_data/文件夹），以及包含 Python 脚本的simulate.py文件，我们可以在命令行运行该脚本来模拟本章节的数据。

模拟登录尝试

由于我们很难从泄露数据中找到登录尝试数据（通常因为其敏感性不予共享），我们将进行模拟。模拟需要对统计建模有深入的理解，能够估计某些事件的概率，并根据需要识别适当的假设进行简化。为了运行模拟，我们将构建一个 Python 包（login_attempt_simulator）来模拟登录过程，要求正确的用户名和密码（没有任何额外的身份验证措施，如双重身份验证），以及一个可以在命令行运行的脚本（simulate.py），我们将在本节中讨论这两个内容。

假设

在我们进入处理模拟的代码之前，需要理解一些假设。进行模拟时不可能控制所有可能的变量，因此我们必须确定一些简化假设以便开始。

模拟器对网站有效用户做出以下假设：

• 有效用户按照泊松过程到达，小时的到达率取决于星期几和一天中的时间。泊松过程将每单位时间内的到达（我们的模拟将使用小时）建模为一个均值为λ（lambda）的泊松分布。到达时间服从指数分布，均值为 1/λ。

• 有效用户从 1 到 3 个 IP 地址（每个连接到互联网的设备的唯一标识符）连接，这些 IP 地址由 4 个随机整数组成，范围为[0, 255]，并由句点分隔。尽管极为不可能，但两个有效用户可能会共享一个 IP 地址。

• 有效用户在输入凭证时不太可能犯很多错误。

  重要提示

  到达时间具有无记忆性特性，这意味着两个连续到达之间的时间不会影响随后的到达时间。

模拟器对黑客做出以下假设：

• 黑客试图避免账户锁定，只测试少量的用户名和密码组合，而不是进行全面的字典攻击（对于每个用户，尝试黑客在字典中所有可能的密码）。然而，他们在尝试之间不会添加延迟。

• 由于黑客不想造成拒绝服务攻击，他们限制攻击的频率，每次只进行一次尝试。

• 黑客知道系统中存在的账户数量，并且对用户名的格式有一定了解，但他们只能猜测具体的用户名。他们会选择尝试猜测所有 133 个用户名，或者其中的某些子集。

• 每次攻击都是独立的，也就是说每次攻击都是由单个黑客执行的，而且一个黑客从不进行多次攻击。

• 黑客不会共享哪些用户名-密码组合是正确的。

• 攻击是随机发生的。

• 每个黑客将使用一个唯一的 IP 地址，该地址与有效用户的 IP 地址生成方式相同。然而，我们的模拟器能够改变这个 IP 地址，这是我们将在第十一章《机器学习异常检测》中探讨的功能，目的是让这个场景变得更具挑战性。

• 虽然可能性极小，但黑客的 IP 地址可能与有效用户相同，甚至黑客可能是有效用户。

我们也在抽象化密码猜测的一些复杂性；相反，我们使用随机数字来决定密码是否被猜测正确——这意味着我们没有考虑网站如何存储密码，可能是明文（希望不是）、哈希值（对明文密码进行不可逆转换，使得无需存储实际密码即可进行验证）或加盐哈希值（有关这方面的文章请参见进一步阅读部分）。实际上，黑客可能会获取存储的密码并离线破解它们（请参阅进一步阅读部分结尾处关于彩虹表的文章），在这种情况下，本章讨论的技术可能不会很有帮助，因为日志中不会记录他们的尝试。请记住，本次模拟中的黑客非常显眼。

login_attempt_simulator 包

这个包比上一章中的stock_analysis包要轻量得多；我们只有三个文件：

login_attempt_simulator
|-- __init__.py
|-- login_attempt_simulator.py
`-- utils.py

接下来，我们将在以下章节中逐一讲解这些文件。请注意，部分文档字符串已被删除以简洁起见；请查看文件本身以获取完整的文档。

辅助函数

让我们从utils.py函数开始讨论，它们是我们模拟器类的辅助工具。首先，我们为模块创建文档字符串，并处理导入：

"""Utility functions for the login attempt simulator."""
import ipaddress
import itertools
import json
import random
import string

接下来，我们定义make_user_base()函数，它为我们的 Web 应用程序创建用户库。该函数通过将英语字母表中的一个小写字母与函数内列表中的每个姓氏结合，来创建一个包含用户名的文件，并添加一些管理员账户；这样就形成了一个包含 133 个账户的用户库。通过写入文件，我们确保每次运行模拟时不需要重新生成，而是可以简单地从中读取数据以便将来进行模拟：

def make_user_base(out_file):
    """Generate a user base and save it to a file."""
    with open(out_file, 'w') as user_base:
        for first, last in itertools.product(
            string.ascii_lowercase, 
            ['smith', 'jones', 'kim', 'lopez', 'brown']
        ): # makes 130 accounts
            user_base.write(first + last + '\n')
        # adds 3 more accounts
        for account in ['admin', 'master', 'dba']: 
            user_base.write(account + '\n')

由于我们需要在模拟器中使用这个用户库，我们还写了一个函数来读取用户库文件并将其转化为列表。get_valid_users() 函数将由 make_user_base() 函数写入的文件重新读取到 Python 列表中：

def get_valid_users(user_base_file):
    """Read in users from the user base file."""
    with open(user_base_file, 'r') as file:
        return [user.strip() for user in file.readlines()]

random_ip_generator() 函数生成随机的 IP 地址，格式为 xxx.xxx.xxx.xxx，其中 x 是 [0, 255] 范围内的整数。我们使用 Python 标准库中的 ipaddress 模块（docs.python.org/3/library/ipaddress.html）来避免分配私有 IP 地址：

def random_ip_generator():
    """Randomly generate a fake IP address."""
    try:
        ip_address = ipaddress.IPv4Address('%d.%d.%d.%d' %
            tuple(random.randint(0, 255) for _ in range(4))
        )
    except ipaddress.AddressValueError:
        ip_address = random_ip_generator()
    return str(ip_address) if ip_address.is_global \
        else random_ip_generator()

每个用户将有几个尝试登录的 IP 地址。assign_ip_addresses() 函数为每个用户映射 1 到 3 个随机 IP 地址，创建一个字典：

def assign_ip_addresses(user_list):
    """Assign users 1-3 fake IP addresses."""
    return {
        user: [
            random_ip_generator()
            for _ in range(random.randint(1, 3))
        ] for user in user_list
    }

save_user_ips() 和 read_user_ips() 函数分别将用户-IP 地址映射保存到 JSON 文件中，并将其重新读取到字典文件中：

def save_user_ips(user_ip_dict, file):
    """Save the user-IP address mapping to a JSON file."""
    with open(file, 'w') as file:
        json.dump(user_ip_dict, file)
def read_user_ips(file):
    """Read in the JSON file of the user-IP address mapping."""
    with open(file, 'r') as file:
        return json.loads(file.read())

提示

Python 标准库有许多有用的模块，尽管我们可能不常用，但它们绝对值得了解。在这里，我们使用 json 模块将字典保存到 JSON 文件，并稍后读取它们。我们使用 ipaddress 模块处理 IP 地址，使用 string 模块获取字母表中的字符，而不需要一一输入。

LoginAttemptSimulator 类

login_attempt_simulator.py 文件中的 LoginAttemptSimulator 类负责执行模拟的重负荷工作，包含所有随机数生成逻辑。像往常一样，我们从模块的文档字符串和导入语句开始：

"""Simulator of login attempts from valid users and hackers."""
import calendar
import datetime as dt
from functools import partial
import math
import random
import string
import numpy as np
import pandas as pd
from .utils import random_ip_generator, read_user_ips

接下来，我们开始定义 LoginAttemptSimulator 类及其文档字符串，同时为存储常量定义一些类变量。我们这样做是为了避免使用魔法数字（即代码中看似没有意义的数字）以及避免在多个地方使用字符串时出现拼写错误。请注意，这些信息仅用于我们的日志；Web 应用程序不会向最终用户显示认证尝试失败的原因（也不应该显示）：

class LoginAttemptSimulator:
    """Simulate login attempts from valid users + attackers."""
    ATTEMPTS_BEFORE_LOCKOUT = 3
    ACCOUNT_LOCKED = 'error_account_locked'
    WRONG_USERNAME = 'error_wrong_username'
    WRONG_PASSWORD = 'error_wrong_password'

重要提示

请注意我们如何使用类变量来存储常量，例如错误信息，这样就不会在代码中犯拼写错误。这样，每次使用这些错误信息时，文本都会保持一致，从而保持数据的整洁。在 Python 中，常量通常采用全大写字母形式（www.python.org/dev/peps/pep-0008/#constants）。

__init__() 方法将处理模拟器的设置，例如从指定的文件读取用户库、初始化日志、存储成功概率，并根据需要确定模拟的开始和结束日期：

    def __init__(self, user_base_json_file, start, end=None, *,
                 attacker_success_probs=[.25, .45],
                 valid_user_success_probs=[.87, .93, .95],
                 seed=None):
        # user, ip address dictionary
        self.user_base = read_user_ips(user_base_json_file)
        self.users = [user for user in self.user_base.keys()]
        self.start = start
        self.end = end if end else self.start + \
            dt.timedelta(days=random.uniform(1, 50))
        self.hacker_success_likelihoods = \
            attacker_success_probs
        self.valid_user_success_likelihoods = \
            valid_user_success_probs
        self.log = pd.DataFrame(columns=[
            'datetime', 'source_ip', 'username',
            'success', 'failure_reason'
        ])
        self.hack_log = \
            pd.DataFrame(columns=['start', 'end', 'source_ip'])
        self.locked_accounts = []
        # set seeds for random numbers from random and numpy:
        random.seed(seed)
        np.random.seed(seed)

_record() 方法将每次尝试的结果追加到日志中，记录其来源的 IP 地址、用户名、时间、是否成功以及失败的原因（如果有的话）：

    def _record(self, when, source_ip, username, success, 
                failure_reason):
        """
        Record the outcome of a login attempt.
        Parameters:
            - when: The datetime of the event.
            - source_ip: IP address the attempt came from.
            - username: The username used in the attempt.
            - success: Whether the attempt succeeded (Boolean).
            - failure_reason: Reason for the failure.
        Returns: 
            None, the `log` attribute is updated.
        """
        self.log = self.log.append({
            'datetime': when, 
            'source_ip': source_ip, 
            'username': username, 
            'success': success, 
            'failure_reason': failure_reason
        }, ignore_index=True)

_attempt_login() 方法处理判断登录尝试是否成功的逻辑：

图 8.1 – 模拟逻辑

我们提供输入正确用户名的概率（username_accuracy）以及每次尝试成功输入密码的概率（success_likelihoods）。尝试次数是允许的最大尝试次数与成功概率列表长度（success_likelihoods）中的最小值。每次尝试的结果会通过 functools 传递给 _record()，它允许我们创建函数，将某些参数固定为特定值（这样就不必不断传递相同的值）：

    def _attempt_login(self, when, source_ip, username,
                       username_accuracy, success_likelihoods):
        """
        Simulates a login attempt, allowing for account
        lockouts, and recording the results.
        Parameters:
            - when: The datetime to start trying.
            - source_ip: IP address the attempt came from. 
            - username: The username being used in the attempt.
            - username_accuracy: Prob. username is correct.
            - success_likelihoods: List of probabilities that 
              password is correct (one per attempt).
        Returns:
            The datetime after trying.
        """
        current = when
        recorder = partial(self._record, source_ip=source_ip)
        if random.random() > username_accuracy:
            correct_username = username
            username = self._distort_username(username)
        if username not in self.locked_accounts:
            tries = len(success_likelihoods)
            for i in range(
                min(tries, self.ATTEMPTS_BEFORE_LOCKOUT)
            ):
                current += dt.timedelta(seconds=1)
                if username not in self.users:
                    recorder(
                        when=current, username=username, 
                        success=False,
                        failure_reason=self.WRONG_USERNAME
                    )
                    if random.random() <= username_accuracy:
                        username = correct_username
                    continue
                if random.random() <= success_likelihoods[i]:
                    recorder(
                        when=current, username=username,
                        success=True, failure_reason=None
                    )
                    break
                else:
                    recorder(
                        when=current, username=username, 
                        success=False,
                        failure_reason=self.WRONG_PASSWORD
                    )
            else:
                if tries >= self.ATTEMPTS_BEFORE_LOCKOUT \
                and username in self.users:
                    self.locked_accounts.append(username)
        else:
            recorder(
                when=current, username=username, success=False,
                failure_reason=self.ACCOUNT_LOCKED
            )
            if random.random() >= .5: # unlock account randomly
                self.locked_accounts.remove(username)
        return current

_valid_user_attempts_login() 和 _hacker_attempts_login() 方法是围绕 _attempt_login() 方法的封装，分别处理有效用户和黑客的概率调整。注意，尽管两者都使用高斯（正态）分布来确定用户名的准确性，但有效用户的分布具有更高的均值和更低的标准差，这意味着他们在尝试登录时更有可能提供正确的用户名。这是因为，虽然有效用户可能会打错字（偶尔发生），但黑客则是在猜测：

    def _hacker_attempts_login(self, when, source_ip,
                               username):
        """Simulates a login attempt from an attacker."""
        return self._attempt_login(
            when=when, source_ip=source_ip, username=username,
            username_accuracy=random.gauss(mu=0.35, sigma=0.5),
            success_likelihoods=self.hacker_success_likelihoods
        )
    def _valid_user_attempts_login(self, when, username):
        """Simulates a login attempt from a valid user."""
        return self._attempt_login(
            when=when, username=username,
            source_ip=random.choice(self.user_base[username]),
            username_accuracy=\
                random.gauss(mu=1.01, sigma=0.01),
            success_likelihoods=\
                self.valid_user_success_likelihoods
        )

当模拟器判断用户名无法正确提供时，它会调用 _distort_username() 方法，该方法会随机决定从有效用户名中省略一个字母或将一个字母替换为另一个字母。虽然黑客因猜测而输入错误的用户名（而不是由于打字错误），我们在这里抽象化处理这一细节，以便使用一个统一的函数来为有效用户和黑客引入用户名错误：

    @staticmethod
    def _distort_username(username):
        """
        Alters the username to allow for wrong username login 
        failures. Randomly removes a letter or replaces a 
        letter in a valid username.
        """
        username = list(username)
        change_index = random.randint(0, len(username) - 1)
        if random.random() < .5: # remove random letter
            username.pop(change_index)
        else: # randomly replace a single letter
            username[change_index] = \
                random.choice(string.ascii_lowercase)
        return ''.join(username)

我们使用 _valid_user_arrivals() 方法来生成给定小时内到达的用户数量和使用泊松分布和指数分布生成的到达间隔时间：

    @staticmethod
    def _valid_user_arrivals(when):
        """
        Static method for simulating Poisson process of 
        arrivals (users wanting to log in). Lambda for the 
        Poisson varies depending upon the day and time of week.
        """
        is_weekday = when.weekday() not in (
            calendar.SATURDAY, calendar.SUNDAY
        )
        late_night = when.hour < 5 or when.hour >= 11
        work_time = is_weekday \
                    and (when.hour >= 9 or when.hour <= 17)
        if work_time:
# hours 9-5 on work days get higher lambda 
            poisson_lambda = random.triangular(1.5, 5, 2.75)
        elif late_night:
            # hours in middle of night get lower lambda
            poisson_lambda = random.uniform(0.0, 5.0)
        else:
            poisson_lambda = random.uniform(1.5, 4.25)
        hourly_arrivals = np.random.poisson(poisson_lambda)
        interarrival_times = np.random.exponential(
            1/poisson_lambda, size=hourly_arrivals
        )
        return hourly_arrivals, interarrival_times

重要提示

我们使用 numpy 而不是 random 来从指数分布中生成随机数，因为我们可以一次请求多个值（每个值对应泊松过程确定的每小时到达数）。此外，random 不提供泊松分布，因此我们需要 numpy。

我们的模拟使用了多种不同的分布，因此查看它们的样子可能会有所帮助。以下子图展示了我们使用的每种分布的示例。注意，泊松分布的绘制方式不同。这是因为泊松分布是离散的。因此，我们通常用它来模拟到达情况——在这里，我们用它来模拟尝试登录的用户到达情况。离散分布有一个概率质量函数（PMF），而不是概率密度函数（PDF）：

图 8.2 – 模拟中使用的分布

_hack() 方法为黑客生成一个随机的 IP 地址，并对给定的用户列表进行暴力攻击：

    def _hack(self, when, user_list, vary_ips):
        """
        Simulate an attack by a random hacker.
        Parameters:
            - when: The datetime to start the attack.
            - user_list: The list of users to try to hack.
            - vary_ips: Whether or not to vary the IP address.
        Returns:
            Initial IP address and the end time for recording.
        """
        hacker_ip = random_ip_generator()
        random.shuffle(user_list)
        for user in user_list:
            when = self._hacker_attempts_login(
                when=when, username=user,
                source_ip=random_ip_generator() if vary_ips \
                    else hacker_ip
            )
        return hacker_ip, when

现在我们已经具备了执行仿真主要部分的功能，我们编写 simulate() 方法将所有内容整合在一起：

    def simulate(self, *, attack_prob, try_all_users_prob,
                 vary_ips):
        """
        Simulate login attempts.
        Parameters:
            - attack_probs: Probability of attack in given hour
            - try_all_users_prob: Prob. hacker will try to 
              guess credentials for all users vs random subset.
            - vary_ips: Whether to vary the IP address.
        """
        hours_in_date_range = math.floor(
            (self.end - self.start).total_seconds() / 60 / 60
        )
        for offset in range(hours_in_date_range + 1):
            current = self.start + dt.timedelta(hours=offset)
            # simulate hacker
            if random.random() < attack_prob:
                attack_start = current \
                    + dt.timedelta(hours=random.random())
                source_ip, end_time = self._hack(
                    when=attack_start,
                    user_list=self.users if \
                        random.random() < try_all_users_prob \
                        else random.sample(
self.users, 
                            random.randint(0, len(self.users))
                    ),
                    vary_ips=vary_ips
                )
                self.hack_log = self.hack_log.append(
                    dict(
                        start=attack_start, end=end_time, 
                        source_ip=source_ip
                    ), ignore_index=True
                )
            # simulate valid users
            hourly_arrivals, interarrival_times = \
                self._valid_user_arrivals(current)
            random_user = random.choice(self.users)
            random_ip = \
                random.choice(self.user_base[random_user])
            for i in range(hourly_arrivals):
                current += \
                    dt.timedelta(hours=interarrival_times[i])
                current = self._valid_user_attempts_login(
                    current, random_user
                )

我们希望将日志保存为 CSV 文件，因此我们将 _save() 方法作为静态方法添加，以减少两个保存方法中的代码重复。save_log() 方法将保存登录尝试，而 save_hack_log() 方法将保存攻击记录：

    @staticmethod
    def _save(data, filename, sort_column):
        """Sort data by the datetime and save to a CSV file."""
        data.sort_values(sort_column)\
            .to_csv(filename, index=False)
    def save_log(self, filename):
        """Save the login attempts log to a CSV file."""
        self._save(self.log, filename, 'datetime')
    def save_hack_log(self, filename):
        """Save the record of the attacks to a CSV file."""
        self._save(self.hack_log, filename, 'start')

注意到这个类中有很多私有方法，这是因为该类的用户只需要能够创建该类的实例（__init__()），按小时进行仿真（simulate()），并保存输出（save_log() 和 save_hack_log()）——所有其他方法仅供该类的对象内部使用。后台的这些方法将处理大部分工作。

最后，我们有 __init__.py 文件，这使得这个目录变成一个包，同时也为我们提供了一种更容易的方式来导入主类：

"""Package for simulating login data."""
from .login_attempt_simulator import LoginAttemptSimulator

现在我们已经了解了仿真器是如何工作的，接下来我们将讨论如何运行仿真以收集登录尝试的数据。

从命令行进行仿真

我们可以将模拟登录尝试的代码包装成一个可以轻松通过命令行运行的脚本，而不是每次都写代码。Python 标准库中有 argparse 模块（docs.python.org/3/library/argparse.html），允许我们为脚本指定可以从命令行提供的参数。

让我们来看看 simulate.py 文件，看看如何做到这一点。我们从导入开始：

"""Script for simulating login attempts."""
import argparse
import datetime as dt
import os
import logging
import random
import login_attempt_simulator as sim

为了在命令行中使用时提供状态更新，我们将使用标准库中的 logging 模块设置日志消息（docs.python.org/3/library/logging.html）：

# Logging configuration
FORMAT = '[%(levelname)s] [ %(name)s ] %(message)s'
logging.basicConfig(level=logging.INFO, format=FORMAT)
logger = logging.getLogger(os.path.basename(__file__))

接下来，我们定义了一些实用函数，用于生成我们在仿真过程中读取和写入数据时需要的文件路径：

def get_simulation_file_path(path_provided, directory,
                             default_file):
    """Get filepath, make directory if necessary."""
    if path_provided:
        file = path_provided
    else:
        if not os.path.exists(directory):
            os.mkdir(directory)
        file = os.path.join(directory, default_file)
    return file
def get_user_base_file_path(path_provided, default_file):
    """Get the path for a user_data directory file."""
    return get_simulation_file_path(
        path_provided, 'user_data', default_file
    )
def get_log_file_path(path_provided, default_file):
    """Get the path for a logs directory file."""
    return get_simulation_file_path(
        path_provided, 'logs', default_file
    )

这个脚本的最大部分定义了可以传递的命令行参数——我们将允许用户指定是否要创建一个新的用户基础，设置种子，仿真开始时间，仿真时长，以及保存所有文件的位置。实际的仿真通过我们构建的包在几行代码内完成。这个部分只有在运行该模块时才会执行，而不是在导入时：

if __name__ == '__main__':
    # command-line argument parsing
    parser = argparse.ArgumentParser()
    parser.add_argument(
        'days', type=float,
        help='number of days to simulate from start'
    )
    parser.add_argument(
        'start_date', type=str,
        help="datetime to start in the form 'YYYY-MM-DD(...)'"
    )
    parser.add_argument(
'-m', '--make', action='store_true', 
        help='make user base'
    )
    parser.add_argument(
'-s', '--seed', type=int, 
        help='set a seed for reproducibility'
    )
    parser.add_argument(
'-u', '--userbase', 
        help='file to write the user base to'
    )
    parser.add_argument(
'-i', '--ip', 
        help='file to write user-IP address map to'
    )
    parser.add_argument(
        '-l', '--log', help='file to write the attempt log to'
    )
    parser.add_argument(
'-hl', '--hacklog', 
        help='file to write the hack log to'
    )

提示

if __name__ == '__main__' 块中放置的代码只有在该模块作为脚本运行时才会被执行。这样，我们就可以在不运行仿真的情况下导入模块中定义的函数。

定义好参数后，我们需要解析它们才能使用：

    args = parser.parse_args()

一旦我们解析了命令行参数，就检查是否需要生成用户基础数据或读取它：

    user_ip_mapping_file = \
        get_user_base_file_path(args.ip, 'user_ips.json')
    if args.make:
        logger.warning(
            'Creating new user base, mapping IP addresses.'
        )
        user_base_file = get_user_base_file_path(
            args.userbase, 'user_base.txt'
        )
        # seed the creation of user base
        random.seed(args.seed)
        # create usernames and write to file
        sim.utils.make_user_base(user_base_file)
        # create 1 or more IP addresses per user, save mapping 
        valid_users = sim.utils.get_valid_users(user_base_file)
        sim.utils.save_user_ips(
            sim.utils.assign_ip_addresses(valid_users), 
            user_ip_mapping_file
        )

之后，我们从命令行参数中解析起始日期，并通过将持续时间添加到起始日期来确定结束日期：

    try:
        start = \
            dt.datetime(*map(int, args.start_date.split('-')))
    except TypeError:
        logger.error('Start date must be in "YYYY-MM-DD" form')
        raise
    except ValueError:
        logger.warning(
            f'Could not interpret {args.start_date}, '
            'using January 1, 2020 at 12AM as start instead'
        )
        start = dt.datetime(2020, 1, 1)
    end = start + dt.timedelta(days=args.days)

提示

try clause and multiple except clauses. We can specify how to handle specific errors occurring during code execution (called except clause. In this case, we have the logger object print a more helpful message for the user, and then re-raise the same exception (because we don't intend to handle it) by simply writing raise. This ends the program—the user can then try again with valid input. Try triggering this exception to see how much more useful this is. One thing to keep in mind, though, is that order matters—be sure to handle specific exceptions before having a general except clause; otherwise, the code specific to each exception type will never trigger. Also, note that using except without providing a specific exception will catch everything, even exceptions not meant to be caught.

最后，我们运行实际的模拟并将结果写入指定的文件（或默认路径）。我们将某一小时内发生攻击的概率设置为 10%（attack_prob），黑客尝试猜测所有用户名的概率设置为 20%（try_all_users_prob），并让黑客在所有尝试中使用相同的 IP 地址（vary_ips）：

    try:
        logger.info(f'Simulating {args.days} days...')
        simulator = sim.LoginAttemptSimulator(
            user_ip_mapping_file, start, end, seed=args.seed
        )
        simulator.simulate(
            attack_prob=0.1, try_all_users_prob=0.2, 
            vary_ips=False
        )
        # save logs
        logger.info('Saving logs')
        simulator.save_hack_log(
            get_log_file_path(args.hacklog, 'attacks.csv')
        )
        simulator.save_log(
            get_log_file_path(args.log, 'log.csv')
        )
        logger.info('All done!')
    except:
        logger.error('Oops! Something went wrong...')
        raise

提示

请注意，我们使用了logger对象在脚本中打印有用的信息到屏幕；这将帮助脚本用户了解进程进行到什么阶段。这些信息有不同的严重级别（我们在这里使用的是INFO、WARNING和ERROR），允许它们用于调试（DEBUG级别），并且在代码进入生产阶段后，打印的最小级别可以提高到INFO，这样就不会再打印DEBUG信息。这比简单的print()语句要强大得多，因为我们不需要担心在进入生产环境时删除它们，或者在开发继续进行时再将这些信息加回去。

现在让我们来看看如何运行这个脚本。我们知道simulate.py可以在命令行上运行，但我们如何查看需要传递哪些参数呢？很简单——我们只需在调用时添加帮助标志（-h或--help）：

(book_env) $ python3 simulate.py -h
usage: simulate.py [-h] [-m] [-s SEED] [-u USERBASE] [-i IP] 
                   [-l LOG] [-hl HACKLOG]
                   days start_date
positional arguments:
  days                  number of days to simulate from start
  start_date            datetime to start in the form 
                        'YYYY-MM-DD' or 'YYYY-MM-DD-HH'
optional arguments:
  -h, --help            show this help message and exit
  -m, --make            make user base
  -s SEED, --seed SEED  set a seed for reproducibility
  -u USERBASE, --userbase USERBASE
                        file to write the user base to
  -i IP, --ip IP        file to write the user-IP address 
                        map to
  -l LOG, --log LOG     file to write the attempt log to
  -hl HACKLOG, --hacklog HACKLOG
                        file to write the hack log to

重要提示

请注意，当我们使用argparse添加其他参数时，并没有指定help参数；它是由argparse自动创建的。

一旦我们知道可以传递哪些参数，并决定了我们想要提供哪些参数，我们就可以运行模拟了。让我们模拟从 2018 年 11 月 1 日 12 点开始的 30 天，并让脚本创建所需的用户基础和 IP 地址映射：

(book_env) $ python3 simulate.py -ms 0 30 '2018-11-01'
[WARNING] [ simulate.py ] Creating new user base and mapping IP addresses to them.
[INFO] [ simulate.py ] Simulating 30.0 days...
[INFO] [ simulate.py ] Saving logs
[INFO] [ simulate.py ] All done!

提示

由于我们设置了种子（-s 0），因此此模拟的输出是可重现的。只需移除种子或更改它，就可以得到不同的结果。

Python 模块也可以作为脚本运行。与导入模块不同，当我们将模块作为脚本运行时，if __name__ == '__main__'下方的任何代码也会被执行，这意味着我们不总是需要编写单独的脚本。我们构建的大多数模块只定义了函数和类，因此作为脚本运行并不会有任何效果；然而，我们在第一章《数据分析简介》中创建虚拟环境时，正是这种方式。前面的代码块因此等同于以下命令：

# leave off the .py
(book_env) $ python3 -m simulate -ms 0 30 "2018-11-01"

现在我们已经有了模拟数据，开始进行分析吧。

探索性数据分析

在这种情况下，我们有幸可以访问标记数据（logs/attacks.csv），并将其用于研究如何区分合法用户和攻击者。然而，这种情况在我们离开研究阶段并进入应用阶段时往往无法实现。在第十一章《机器学习异常检测》中，我们将重新审视这一场景，但从没有标记数据开始，以增加挑战。像往常一样，我们先进行导入并读取数据：

>>> %matplotlib inline
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> import pandas as pd
>>> import seaborn as sns
>>> log = pd.read_csv(
...     'logs/log.csv', index_col='datetime', parse_dates=True
... )

登录尝试数据框（log）包含了每次尝试的日期和时间（datetime列）、来源的 IP 地址（source_ip）、使用的用户名（username）、是否成功（success）以及如果失败的话失败原因（failure_reason）：

图 8.3 – 登录尝试数据示例

在处理这些数据时，我们需要思考什么是正常活动，什么是黑客活动。两者之间的任何显著差异都可能被用来识别黑客。我们预计合法用户的成功率较高，最常见的失败原因是密码错误。我们预计用户会从不同的 IP 地址登录（手机、家用电脑、工作电脑以及其他可能的设备），并且有可能人们共享设备。由于我们不了解该网络应用程序的性质，我们不能说多次登录是否正常。我们也不知道这些数据的时区，因此不能对登录时间做出任何推断。理论上，我们可以查看这些 IP 地址来自哪些国家，但由于有方法可以隐藏 IP 地址，因此我们不会走这条路。基于现有数据，我们可以选择以下几种可行的方式：

• 调查登录尝试和失败的任何异常波动（无论是总体还是按 IP 地址统计）。

• 检查失败原因是用户名错误的情况。

• 查看每个 IP 地址的失败率。

• 查找尝试使用多个不同用户名登录的 IP 地址。

另一个需要注意的点是，我们希望尽早标记异常行为，而不是等到最后。等待一个月才标记某个行为的价值较低（随着时间推移，价值迅速下降），因此我们需要找到更早标记异常的方式；例如，使用每小时的频率。由于我们处于研究阶段，可以使用一些标记数据：

>>> attacks = pd.read_csv(
...     'logs/attacks.csv',
...     converters={
...         'start': np.datetime64, 
...         'end': np.datetime64
...     }
... ) # make start and end columns datetimes but not the index

这些数据是针对网络应用程序（attacks）的攻击记录。它包含了攻击开始的日期和时间（start）、攻击结束的日期和时间（end）以及与攻击相关的 IP 地址（source_ip）：

图 8.4 – 标记数据示例

使用shape属性，我们可以看到 72 次攻击和 12,836 次来自有效用户和恶意用户的登录尝试，使用nunique()，我们看到 22%的 IP 地址与攻击相关：

>>> attacks.shape, log.shape
((72, 3), (12836, 4))
>>> attacks.source_ip.nunique() / log.source_ip.nunique()
0.22018348623853212

重要提示

通常情况下，仅凭这些数据很难知道攻击发生的具体时间——攻击可以在不被发现的情况下持续很长时间，即使是这样，也并不容易将攻击者的行为与正常用户的行为区分开来。

我们的数据相当干净（毕竟我们就是为这个目的设计的），所以让我们看看通过执行一些探索性数据分析（EDA）是否能发现有趣的东西。首先，让我们看看每小时的登录尝试数量：

>>> log.assign(attempts=1).attempts.resample('1H').sum()\
...     .plot(figsize=(15, 5), title='hourly attempts')\
...     .set(xlabel='datetime', ylabel='attempts')

几个小时的登录尝试出现了非常大的峰值，这可能是攻击发生的时间。使用这张图，我们可以报告登录尝试活动较高的小时数，但仅此而已：

图 8.5 – 每小时登录尝试

另一个有趣的探索方向是查看每个 IP 地址的尝试次数。我们可以通过运行以下命令来实现这一点：

>>> log.source_ip.value_counts().describe()
count    327.000000
mean      39.253823
std       69.279330
min        1.000000
25%        5.000000
50%       10.000000
75%       22.500000
max      257.000000
Name: source_ip, dtype: float64

这些数据显然有一些异常值，它们将每个 IP 地址的尝试次数拉得很高。让我们创建一些图表来更好地评估这一点：

>>> fig, axes = plt.subplots(1, 2, figsize=(15, 5))
>>> log.source_ip.value_counts()\
...     .plot(kind='box', ax=axes[0]).set_ylabel('attempts')
>>> log.source_ip.value_counts()\
...     .plot(kind='hist', bins=50, ax=axes[1])\
...     .set_xlabel('attempts')
>>> fig.suptitle('Attempts per IP Address')

每个 IP 地址的尝试分布是有效用户和攻击者分布的总和。直方图显示这个分布是双峰的，但仅凭这张图我们无法确定所有高尝试次数的 IP 地址是否都是黑客：

图 8.6 – 每个 IP 地址的登录尝试分布

由于我们可以访问每次攻击的详细信息，我们可以检查直方图的右侧部分是否是黑客的分布。根据尝试次数，黑客的 IP 地址占前列 IP 地址的 88.9%：

>>> num_hackers = attacks.source_ip.nunique()
>>> log.source_ip.value_counts().index[:num_hackers]\
...     .isin(attacks.source_ip).sum() / num_hackers
0.8888888888888888

我们可以在此停止，并标记出每个月尝试次数最多的 IP 地址列表，但我们更可能需要一个更具鲁棒性的解决方案，因为黑客每次都可能更改 IP 地址来避免检测。理想情况下，我们还希望在不等待完整一个月的数据的情况下就能检测到攻击。然而，通过查看每个 IP 地址的每小时尝试次数，遗憾的是，似乎并没有提供太多信息：

>>> log.assign(attempts=1).groupby('source_ip').attempts\
...     .resample('1H').sum().unstack().mean()\
...     .plot(
...         figsize=(15, 5), 
...         title='average hourly attempts per IP address'
...     ).set_ylabel('average hourly attempts per IP address')

请回想一下来自第一章的内容，数据分析简介，其中提到均值对异常值不具鲁棒性。如果攻击者进行了多次尝试，异常值会将每个 IP 地址的平均每小时登录尝试次数拉高。我们可以在这条线图中看到几个大的峰值，但注意到它们中的许多只有两到三次。我们真的能期望一个用户只通过一个 IP 地址访问 Web 应用程序吗？这很可能不是一个现实的假设：

图 8.7 – 每个 IP 地址的平均每小时登录尝试次数

所以，如果我们不能仅依赖 IP 地址（毕竟，黑客可能足够聪明，能够将攻击分布到多个不同的地址上），我们还能尝试什么呢？或许黑客在成功登录方面遇到了更多困难：

>>> log[log.source_ip.isin(attacks.source_ip)]\
...     .success.value_counts(normalize=True)
False    0.831801
True     0.168199
Name: success, dtype: float64

黑客成功的概率只有 17%，但合法用户的成功率有多高呢？这些信息对于确定网站正常行为的基准值非常重要。正如我们所预期的，合法用户的成功率要高得多：

>>> log[~log.source_ip.isin(attacks.source_ip)]\
...     .success.value_counts(normalize=True)
True     0.873957
False    0.126043
Name: success, dtype: float64

由于日志中包含了登录失败的原因，我们可以使用交叉表来查看黑客和合法用户未能成功登录的原因。这里的任何差异都可能帮助我们区分这两组用户：

>>> pd.crosstab(
...     index=pd.Series(
...         log.source_ip.isin(attacks.source_ip),
...         name='is_hacker'
...     ), columns=log.failure_reason
... )

合法用户有时会输入错误的密码或用户名，但黑客在同时输入正确的用户名和密码时会遇到更多问题：

图 8.8 – 登录失败尝试的原因

合法用户在输入凭证时不会犯太多错误，因此如果黑客尝试登录的次数很多，并且涉及多个用户，我们可以标记这种行为。为了确认这一点，我们可以查看每个用户的平均每小时尝试次数：

>>> log.assign(attempts=1).groupby('username').attempts\
...     .resample('1H').sum().unstack().mean()\
...     .plot(figsize=(15, 5),
...           title='average hourly attempts per user')\
...     .set_ylabel('average hourly attempts per user')

大多数情况下，每个用户名每小时的尝试次数不到一次。这个指标的波动也不能保证就是攻击的迹象。或许是网站正在进行闪购活动；在这种情况下，我们很可能会看到由合法用户引起的这个指标的激增：

图 8.9 – 每个用户名的平均每小时登录尝试次数

根据我们的发现，错误率似乎是检测攻击的最有用指标，因此我们将研究错误率高的 IP 地址。为此，我们可以创建一个透视表来计算一些有用的指标：

>>> pivot = log.pivot_table(
...     values='success', index=log.source_ip, 
...     columns=log.failure_reason.fillna('success'), 
...     aggfunc='count', fill_value=0
... )
>>> pivot.insert(0, 'attempts', pivot.sum(axis=1))
>>> pivot = pivot.sort_values('attempts', ascending=False)\
...     .assign(
...         success_rate=lambda x: x.success / x.attempts,
...         error_rate=lambda x: 1 - x.success_rate
...     )
>>> pivot.head()

提示

insert()方法允许我们将新创建的attempts列插入到当前数据框的特定位置，并且是就地操作。我们创建了attempts列，它是错误和成功的总和（我们将failure_reason列中的NaN值用success填充，以便在此处进行计数），通过在axis=1方向上进行求和。

这生成了以下按尝试次数排序的透视表（从最多到最少）：

图 8.10 – 每个 IP 地址的度量

我们知道某些 IP 地址正在进行多次尝试，因此值得调查每个 IP 地址尝试登录的用户名数量；我们预计合法用户只会从少数几个 IP 地址登录，并且不会与许多其他用户共享他们的 IP 地址。这可以通过分组和聚合来确定：

>>> log.groupby('source_ip').agg(dict(username='nunique'))\
...     .username.value_counts().describe()
count     53.000000
mean       6.169811
std       34.562505
min        1.000000
25%        1.000000
50%        1.000000
75%        2.000000
max      253.000000
Name: username, dtype: float64

这看起来确实是隔离恶意用户的好策略。大多数 IP 地址由两名或更少的用户使用，但最大值为 253 个。虽然这个标准可能有助于我们识别一些攻击者，但如果黑客足够聪明，可以在攻击过程中不断更换 IP 地址，这个标准就无法起到作用了。

在我们继续讨论异常检测方法之前，先看看我们是否能通过视觉识别出黑客。让我们为每个 IP 地址的成功次数和尝试次数绘制散点图：

>>> pivot.plot(
...     kind='scatter', x='attempts', y='success', alpha=0.25,
...     title='successes vs. attempts by IP address' 
... )

似乎有几个明显的聚类。在图的左下角，我们看到一些点形成了一条成功与尝试呈一对一关系的线。右上部分则包含一个较为稀疏的聚类，尝试次数较高，成功次数适中。由于我们使用了alpha参数来控制透明度，我们可以看到，似乎连接这两个聚类的点迹并不密集。即便没有坐标轴的比例尺，我们也能预测左下角的聚类是普通用户，右上角的聚类是黑客（因为我们假设普通用户比黑客多，而且普通用户的成功率较高）。然而，中间的点则更难判断：

图 8.11 – 每个 IP 地址的成功与尝试的散点图

在不做任何假设的情况下，我们可以绘制一条边界线，将中间的点与其最近的聚类分组：

>>> ax = pivot.plot(
...     kind='scatter', x='attempts', y='success', alpha=0.25, 
...     title='successes vs. attempts by IP address'
... )
>>> plt.axvline(
...     125, label='sample boundary',
...     color='red', linestyle='--'
... )
>>> plt.legend(loc='lower right')

当然，在缺乏标注数据的情况下，评估这个决策边界的有效性是困难的：

图 8.12 – 可视化决策边界

幸运的是，我们有黑客使用的 IP 地址数据，因为我们已经获得了标注数据来进行研究，所以我们可以使用seaborn来实际看到这种分离：

>>> fig, axes = plt.subplots(1, 2, figsize=(15, 5))
>>> for ax in axes:
...     sns.scatterplot(
...         y=pivot.success, x=pivot.attempts, 
...         hue=pivot.assign(
...             is_hacker=\
...                 lambda x: x.index.isin(attacks.source_ip)
...         ).is_hacker,
...         ax=ax, alpha=0.5
...     ) 
...     for spine in ['top', 'right']: # make less boxy
...         ax.spines[spine].set_visible(False)
>>> axes[1].set_xscale('log')
>>> plt.suptitle('successes vs. attempts by IP address')

我们关于存在两个明显聚类的直觉完全正确。然而，中间区域的判断则要复杂得多。左侧的蓝色（较深）点似乎呈向上排列，而左侧的橙色（较浅）点则跟随一条通向橙色聚类的线条。通过绘制尝试次数的对数，我们能更好地将橙色中间点和蓝色点分开：

图 8.13 – 使用标注数据来验证我们的直觉

记住，我们还可以使用箱型图来检查是否存在异常值，异常值会显示为点。让我们看看每个 IP 地址的成功与尝试分布情况：

>>> pivot[['attempts', 'success']].plot(
...     kind='box', subplots=True, figsize=(10, 3),
...     title='stats per IP address'
... )

被标记为离群值的点与我们之前绘制的散点图右上角的点重合：

图 8.14 – 检查异常值

现在我们已经对数据有了充分的了解，我们准备学习如何实现一些简单的异常检测策略。

实现基于规则的异常检测

是时候抓住那些黑客了。在前一节的 EDA（探索性数据分析）之后，我们已经对如何进行此操作有了一定的了解。实际上，这要困难得多，因为它涉及更多的维度，但我们在这里简化了处理过程。我们希望找到那些尝试次数过多但成功率较低的 IP 地址，以及那些使用比我们认为正常的更多独特用户名进行登录尝试的 IP 地址（异常行为）。为了实现这一点，我们将采用基于阈值的规则作为我们进行异常检测的第一步；接下来，在第十一章《机器学习异常检测》中，我们将探讨一些机器学习技术，并重新审视这一场景。

由于我们有兴趣标记可疑的 IP 地址，我们将安排数据，以便每个 IP 地址有每小时的聚合数据（如果该小时内有活动）：

>>> hourly_ip_logs = log.assign(
...     failures=lambda x: np.invert(x.success)
... ).groupby('source_ip').resample('1H').agg({
...     'username': 'nunique', 'success': 'sum', 
...     'failures': 'sum'
... }).assign(
...     attempts=lambda x: x.success + x.failures,
...     success_rate=lambda x: x.success / x.attempts,
...     failure_rate=lambda x: 1 - x.success_rate
... ).dropna().reset_index()

提示

np.invert()函数是一个轻松翻转布尔值的方法。它将True转换为False，将False转换为True，适用于 NumPy 数组结构。

聚合后的数据如下所示：

图 8.15 – 每个 IP 地址的每小时聚合数据

最简单的基于规则的异常检测方法是计算阈值，并检查数据是否超出该阈值。这可能意味着值低于某个下限阈值，或超过某个上限阈值。由于我们关注的是登录尝试，我们对高于正常水平的值感兴趣。因此，我们将计算上限阈值，并将其与我们的数据进行比较。

百分比差异

假设我们对网站上正常的登录尝试活动（排除黑客的影响）有一定了解，我们可以通过某个百分比的偏差来标记那些与此偏离的值。为了计算这个基准，我们可以随机抽取一些 IP 地址（每小时重复抽取），并计算它们的平均登录尝试次数。由于数据量较少（每个小时大约有 50 个独特的 IP 地址可供选择），我们采用自助法（bootstrap）。

为了实现这一点，我们可以编写一个函数，该函数接受我们刚刚创建的聚合数据框，并传入每一列的数据统计名称作为计算阈值的起始点：

>>> def get_baselines(hourly_ip_logs, func, *args, **kwargs):
...     """
...     Calculate hourly bootstrapped statistic per column.
...
...     Parameters:
...         - hourly_ip_logs: Data to sample from.
...         - func: Statistic to calculate.
...         - args: Additional positional arguments for `func`
...         - kwargs: Additional keyword arguments for `func`
...
...     Returns: 
...         `DataFrame` of hourly bootstrapped statistics
...     """
...     if isinstance(func, str):
...         func = getattr(pd.DataFrame, func)
...
...     return hourly_ip_logs.assign(
...         hour=lambda x: x.datetime.dt.hour
...     ).groupby('hour').apply(
...         lambda x: x\
...             .sample(10, random_state=0, replace=True)\
...             .pipe(func, *args, **kwargs, numeric_only=True)
...     )

重要提示

random_state is used with sample() for reproducibility; however, in practice, we will probably not want to always pick the same rows.

请注意，如果我们在apply()内使用sample()，在按我们想要抽样的列分组后，我们可以为所有组（这里是小时）获得大小相等的样本。这意味着我们在每小时对每列进行有放回的抽样，选择 10 行。我们必须按小时进行抽样，因为如果进行简单的随机抽样，可能会没有每小时的统计数据。让我们使用get_baselines()通过均值计算列基准：

>>> averages = get_baselines(hourly_ip_logs, 'mean')
>>> averages.shape
(24, 7)

提示

如果我们想进行分层随机抽样，可以将get_baselines()函数中的10替换为x.shape[0] * pct，其中pct是我们希望从每个组中抽取的百分比。

每列都包含通过随机选择的 10 个 IP 地址估算正常行为的每小时均值。然而，这种方法并不能保证我们不会将黑客活动混入基准计算中。例如，我们来看一下故障率基准值最高的六个小时：

>>> averages.nlargest(6, 'failure_rate')

我们可能会发现很难在19、23或14点钟时标记任何活动为异常，因为这些小时的故障率和尝试过的独特用户名都很高：

图 8.16 – 使用均值计算的每小时基准

为了应对这个问题，我们可以通过让排名前* x *%的值在基准计算中无效来修剪我们的摘要统计信息。我们将从每个小时的数据中移除超过第 95 百分位的数据。首先，我们编写一个函数，用于修剪某个小时内超出给定分位数的行：

>>> def trim(x, quantile):
...     """
...     Remove rows with entries for the username, attempts, 
...     or failure_rate columns above a given quantile.
...     """
...     mask = (
...         (x.username <= x.username.quantile(quantile)) &
...         (x.attempts <= x.attempts.quantile(quantile)) &
...         (x.failure_rate
...          <= x.failure_rate.quantile(quantile))
...     )
...     return x[mask]

接下来，我们将按小时对 IP 地址数据进行分组，并应用我们的修剪函数。由于我们将使用自举函数，因此需要清理一些由此操作产生的多余列，所以我们删除hour列，重置索引，然后移除分组列和旧的索引：

>>> trimmed_hourly_logs = hourly_ip_logs\
...     .assign(hour=lambda x: x.datetime.dt.hour)\
...     .groupby('hour').apply(lambda x: trim(x, 0.95))\
...     .drop(columns='hour').reset_index().iloc[:,2:]

现在，我们可以使用get_baselines()函数，通过平均值和修剪后的数据来获取我们的基准：

>>> averages = get_baselines(trimmed_hourly_logs, 'mean')
>>> averages.iloc[[19, 23, 3, 11, 14, 16]]

经过修剪后的基准在19、23和14小时与图 8.16相比有了显著变化：

图 8.17 – 使用均值修剪后的每小时基准

现在我们已经有了基准，接下来我们来写一个函数，负责从基准和每列的百分比差异中计算阈值，并返回被标记为黑客的 IP 地址：

>>> def pct_change_threshold(hourly_ip_logs, baselines,
...                          pcts=None):
...     """
...     Return flagged IP addresses based on thresholds.
...
...     Parameters:
...         - hourly_ip_logs: Aggregated data per IP address.
...         - baselines: Hourly baselines per column in data.
...         - pcts: Dictionary of custom percentages per column 
...           for calculating upper bound thresholds
...           (baseline * pct). If not provided, pct will be 1
...
...     Returns: `Series` containing the IP addresses flagged.
...     """
...     pcts = {} if not pcts else pcts
...
...     return hourly_ip_logs.assign(
...         hour=lambda x: x.datetime.dt.hour
...     ).join(
...         baselines, on='hour', rsuffix='_baseline'
...     ).assign(
...         too_many_users=lambda x: x.username_baseline \
...             * pcts.get('username', 1) <= x.username,
...         too_many_attempts=lambda x: x.attempts_baseline \
...             * pcts.get('attempts', 1) <= x.attempts,
...         high_failure_rate=lambda x: \
...             x.failure_rate_baseline \
...             * pcts.get('failure_rate', 1) <= x.failure_rate
...     ).query(
...         'too_many_users and too_many_attempts '
...         'and high_failure_rate'
...     ).source_ip.drop_duplicates()

pct_change_threshold()函数使用一系列链式操作来返回被标记的 IP 地址：

1. 首先，它将基准与hour列上的每小时 IP 地址日志连接。由于所有基准列的名称与每小时 IP 地址日志相同，并且我们不希望根据这些名称进行连接，因此我们在列名后添加'_baseline'后缀。

2. 之后，所有需要检查是否超出阈值的数据都在同一个数据框中。我们使用assign()方法创建三个新的布尔列，表示我们每个条件（过多用户、过多尝试和高失败率）是否被违反。

3. 然后，我们链式调用query()方法，这样我们就可以轻松选择所有布尔列值为True的行（注意我们无需显式地写出<column> == True）。

4. 最后，我们确保只返回 IP 地址，并删除任何重复项，以防相同的 IP 地址在多个小时内被标记。

为了使用这个函数，我们需要选择每个基线的百分比差异。默认情况下，这将是基线的 100%，而基线作为平均值，可能会标记出过多的 IP 地址。因此，让我们选择将每个标准的 IP 地址提高 25%的阈值：

>>> pct_from_mean_ips = pct_change_threshold(
...     hourly_ip_logs, averages, 
...     {key: 1.25 for key in [
...         'username', 'attempts', 'failure_rate'
...     ]}
... )

提示

我们使用的百分比存储在字典中，字典的键是对应列名，值是百分比。如果函数调用者没有提供这些值，我们会使用默认的 100%，因为我们使用get()方法从字典中选择。

这些规则标记了 73 个 IP 地址：

>>> pct_from_mean_ips.nunique()
73

重要提示

在实际操作中，我们可能不会对用于计算基线的条目运行此规则，因为它们的行为会影响基线的定义。

Tukey fence

正如我们在第一章《数据分析简介》中讨论的那样，均值对离群值并不稳健。如果我们觉得有很多离群值影响了我们的基线，我们可以返回到百分比差异，尝试使用中位数，或考虑使用Tukey fence。记住，在之前的章节中，我们提到 Tukey fence 的边界来自第一四分位数和第三四分位数，以及四分位间距（IQR）。因为我们只关心超出上限的值，这就解决了均值的问题，前提是离群值占数据量的比例小于 25%。我们可以使用以下公式来计算上限：

我们的get_baselines()函数仍然会帮助我们，但我们需要做一些额外的处理。我们将编写一个函数来计算 Tukey fence 的上限，并让我们测试不同的乘数（k）值。注意，我们在这里也可以选择使用 Tukey fence 的百分比：

>>> def tukey_fence_test(trimmed_data, logs, k, pct=None):
...     """
...     See which IP addresses get flagged with a Tukey fence 
...     with multiplier k and optional percent differences.
...  
...     Parameters: 
...         - trimmed_data: Data for calculating the baselines
...         - logs: The data to test
...         - k: The multiplier for the IQR
...         - pct: Dictionary of percentages per column for use 
...                with `pct_change_threshold()`
...
...     Returns: 
...         `pandas.Series` of flagged IP addresses
...     """
...     q3 = get_baselines(trimmed_data, 'quantile', .75)\
...         .drop(columns=['hour'])
...
...     q1 = get_baselines(trimmed_data, 'quantile', .25)\
...         .drop(columns=['hour'])
...
...     iqr = q3 - q1
...     upper_bound = (q3 + k * iqr).reset_index()
...
...     return pct_change_threshold(logs, upper_bound, pct)

让我们使用tukey_fence_test()函数，利用3的 IQR 乘数，抓取超出 Tukey fence 上限的 IP 地址：

>>> tukey_fence_ips = tukey_fence_test(
...     trimmed_hourly_logs, hourly_ip_logs, k=3
... )

使用这种方法，我们标记了 83 个 IP 地址：

>>> tukey_fence_ips.nunique()
83

重要提示

我们在这里使用了 3 的乘数。然而，根据应用场景，我们可能会看到使用 1.5，以便更加宽松。实际上，我们可以使用任何数字；找到最佳值可能需要一些试错过程。

Z 得分

记住，在第一章《数据分析导论》中，我们还可以计算 Z 分数，并标记距离均值一定标准差的 IP 地址。我们之前写的pct_change_threshold()函数在这种情况下不适用，因为我们不仅仅是在与基线进行比较。相反，我们需要从所有值中减去基线的均值，并除以基线的标准差，因此我们必须重新设计我们的方法。

让我们编写一个新函数z_score_test()，通过任何标准差数量的均值上方作为截断值来执行 Z 分数测试。首先，我们将使用get_baselines()函数，通过修剪后的数据计算每小时的基线标准差。然后，我们将标准差和均值合并，添加后缀。这样，我们就可以将pct_change_threshold()的逻辑应用于这个任务：

>>> def z_score_test(trimmed_data, logs, cutoff):
...     """
...     See which IP addresses get flagged with a Z-score
...     greater than or equal to a cutoff value.
...
...     Parameters: 
...         - trimmed_data: Data for calculating the baselines
...         - logs: The data to test
...         - cutoff: Flag row when z_score >= cutoff
...
...     Returns: 
...         `pandas.Series` of flagged IP addresses
...     """
...     std_dev = get_baselines(trimmed_data, 'std')\
...         .drop(columns=['hour'])
...     averages = get_baselines(trimmed_data, 'mean')\
...         .drop(columns=['hour'])
...
...     return logs.assign(hour=lambda x: x.datetime.dt.hour)\
...         .join(std_dev.join(
...             averages, lsuffix='_std', rsuffix='_mean'
...         ), on='hour')\
...         .assign(
...             too_many_users=lambda x: (
...                 x.username - x.username_mean
...             )/x.username_std >= cutoff,
...             too_many_attempts=lambda x: (
...                 x.attempts - x.attempts_mean
...             )/x.attempts_std >= cutoff,
...             high_failure_rate=lambda x: (
...                 x.failure_rate - x.failure_rate_mean
...             )/x.failure_rate_std >= cutoff
...         ).query(
...             'too_many_users and too_many_attempts '
...             'and high_failure_rate'
...         ).source_ip.drop_duplicates()

让我们使用一个截断值为距离均值三倍标准差或更多的函数：

>>> z_score_ips = \
...     z_score_test(trimmed_hourly_logs, hourly_ip_logs, 3)

使用此方法，我们标记了 62 个 IP 地址：

>>> z_score_ips.nunique()
62

重要提示

在实际应用中，Z 分数的截断值也是我们需要调整的一个参数。

性能评估

因此，我们现在有一系列的 IP 地址，分别对应于每组规则，但我们希望了解每种方法的效果如何（假设我们可以实际检查）。在这种情况下，我们拥有用于研究的攻击者 IP 地址，因此可以查看每种方法标记了多少个正确的地址——这在实践中并不简单；相反，我们可以标记过去发现的恶意地址，并在未来关注类似的行为。

这是一个包含两类的分类问题；我们希望将每个 IP 地址分类为有效用户或恶意用户。这使得我们有四种可能的结果，可以通过混淆矩阵来可视化：

图 8.18 – 混淆矩阵

在此应用中，这些结果的含义如下：

• 真正例 (TP)：我们的方法标记它为恶意的，且它确实是恶意的。

• 真负例 (TN)：我们的方法没有标记它，而且它不是恶意的。

• 假正例 (FP)：我们的方法标记它为恶意的，但它其实不是恶意的。

• 假负例 (FN)：我们的方法没有标记它，但它实际上是恶意的。

真正例和真负例表明我们的方法效果不错，但假正例和假负例是可能需要改进的地方（请记住，这永远不会是完美的）。现在让我们编写一个函数，帮助我们确定每种方法的表现：

>>> def evaluate(alerted_ips, attack_ips, log_ips):
...     """
...     Calculate true positives (TP), false positives (FP),
...     true negatives (TN), and false negatives (FN) for 
...     IP addresses flagged as suspicious.
...
...     Parameters:
...         - alerted_ips: `Series` of flagged IP addresses
...         - attack_ips: `Series` of attacker IP addresses
...         - log_ips: `Series` of all IP addresses seen
...
...     Returns:
...         Tuple of the form (TP, FP, TN, FN)
...     """
...     tp = alerted_ips.isin(attack_ips).sum()
...     tn = np.invert(np.isin(
...         log_ips[~log_ips.isin(alerted_ips)].unique(),
...         attack_ips
...     )).sum()
...     fp = np.invert(alerted_ips.isin(attack_ips)).sum()
...     fn = np.invert(attack_ips.isin(alerted_ips)).sum()
...     return tp, fp, tn, fn

在我们开始计算指标之前，先写一个部分函数，这样我们就不需要不断传递攻击者 IP 地址（attacks.source_ip）和日志中的 IP 地址（pivot.index）。记住，部分函数允许我们固定某些参数的值，然后稍后调用该函数：

>>> from functools import partial
>>> scores = partial(
...     evaluate, attack_ips=attacks.source_ip,
...     log_ips=pivot.index
... )

现在，让我们使用这个来计算一些指标以衡量我们的表现。一个常见的指标是假阳性率（FPR），它告诉我们假警报率。它是通过将假阳性与实际为负的所有情况的比例计算得出的：

假发现率（FDR）是另一种看待假警报的方式，它告诉我们不正确的正例的百分比：

让我们来看看我们使用与均值差异百分比的方法计算的假阳性率（FPR）和假发现率（FDR）：

>>> tp, fp, tn, fn = scores(pct_from_mean_ips)
>>> fp / (fp + tn), fp / (fp + tp)
(0.00392156862745098, 0.0136986301369863)

另一个感兴趣的指标是假阴性率（FNR），它告诉我们我们未能检测到的内容（漏检率）。它通过将假阴性与实际为正的所有情况的比例计算得出：

观察假阴性的一种替代方式是假遗漏率（FOR），它告诉我们错误标记为负例的案例百分比：

我们的均值差异百分比方法没有假阴性，因此 FNR 和 FOR 都是零：

>>> fn / (fn + tp), fn / (fn + tn)
(0.0, 0.0)

这里通常存在一个权衡——我们是否想尽可能捕捉到更多的黑客活动，并冒着标记有效用户的风险（通过关注 FNR/FOR），还是我们希望避免给有效用户带来不便，冒着错过黑客活动的风险（通过最小化 FPR/FDR）？这些问题很难回答，将取决于具体领域，因为假阳性的代价不一定等于（甚至远低于）假阴性的代价。

提示

我们将在第九章中讨论更多可以用于评估我们性能的指标，《Python 中的机器学习入门》。

现在让我们编写一个函数来处理所有这些计算：

>>> def classification_stats(tp, fp, tn, fn):
...     """Calculate metrics"""
...     return {
...         'FPR': fp / (fp + tn), 'FDR': fp / (fp + tp),
...         'FNR': fn / (fn + tp), 'FOR': fn / (fn + tn)
...     }

我们现在可以使用evaluate()函数的结果来计算我们的指标。对于均值差异百分比的方法，我们得到以下输出：

>>> classification_stats(tp, fp, tn, fn)
{'FPR': 0.00392156862745098, 'FDR': 0.0136986301369863,
 'FNR': 0.0, 'FOR': 0.0}

看起来我们的三个标准表现得相当不错。如果我们担心在计算基准时会选择黑客 IP 地址，但又不想修剪数据，我们本可以使用中位数代替均值来运行：

>>> medians = get_baselines(hourly_ip_logs, 'median')
>>> pct_from_median_ips = pct_change_threshold(
...     hourly_ip_logs, medians, 
...     {key: 1.25 for key in
...      ['username', 'attempts', 'failure_rate']}
... )

使用中位数时，我们达到了与均值相似的表现。然而，在这种情况下，我们不需要提前修剪数据。这是因为中位数对离群值具有鲁棒性，意味着选择某个小时内的单个黑客 IP 地址并不会像均值那样影响该小时的基准：

>>> tp, fp, tn, fn = scores(pct_from_median_ips)
>>> classification_stats(tp, fp, tn, fn)
{'FPR': 0.00784313725490196, 'FDR': 0.02702702702702703,
 'FNR': 0.0, 'FOR': 0.0}

为了比较讨论过的每种方法，我们可以使用字典推导式来填充一个DataFrame对象，包含性能指标：

>>> pd.DataFrame({
...     method: classification_stats(*scores(ips))
...     for method, ips in {
...         'means': pct_from_mean_ips,
...         'medians': pct_from_median_ips,
...         'Tukey fence': tukey_fence_ips,
...         'Z-scores': z_score_ips
...     }.items()
... })

提示

scores()函数返回一个元组(tp, fp, tn, fn)，但classification_stats()函数期望四个参数。然而，由于score()以classification_stats()期望的顺序返回这些值，我们可以使用*来解包元组并将这些值作为四个位置参数传递。

均值受异常值的影响，但一旦我们对数据进行修剪，它就成了一个可行的方法。我们不需要修剪数据就能处理中位数；中位数的有效性取决于数据中异常值的比例低于 50%。Tukey 围栏法进一步通过使用第三四分位数，并假设数据中异常值的比例低于 25%，将此方法做得更为严格。Z-score 方法也会受到异常值的影响，因为它使用均值；然而，经过修剪的数据使我们能够通过适当的三倍标准差阈值实现良好的性能：

图 8.19 – 性能比较

最终，我们在实际应用中使用哪种方法将取决于假阳性与假阴性带来的成本——是当什么都没有问题时发出警报更糟，还是当有问题时保持沉默更糟？在这种情况下，我们倾向于尽量减少假阴性，因为我们不想错过任何重要信息。

重要提示

异常检测的另一个常见应用是工业环境中的质量或过程控制，例如监控工厂设备的性能和产量。过程控制使用基于阈值和基于模式的规则来判断系统是否失控。这些方法可以用于判断基础数据的分布是否发生了变化，这可能是后续问题的前兆。西电规则和纳尔逊规则是常见的规则。两者的参考资料可以在本章末的进一步阅读部分找到。

总结

在我们的第二个应用章节中，我们学习了如何在 Python 中模拟事件，并额外接触了编写包的技巧。我们还学习了如何编写可以从命令行运行的 Python 脚本，这些脚本被用来运行我们的登录尝试数据模拟。接着，我们对模拟数据进行了探索性数据分析（EDA），以查看是否能够找出哪些因素使得黑客活动容易被识别。

这让我们聚焦于每小时每个 IP 地址尝试验证的不同用户名数量，以及尝试次数和失败率。利用这些指标，我们能够绘制出一个散点图，图中显示了两个不同的点群，另外还有一些连接两个点群的点；显然，这些点代表了有效用户和恶意用户群体，其中一些黑客并不像其他人那样显眼。

最后，我们开始制定规则，以标记出因可疑活动而被怀疑为黑客的 IP 地址。首先，我们使用pandas将数据按小时汇总成每个 IP 地址的数据。然后，我们编写函数，修剪出大于 95 百分位的数据，并计算每小时给定统计数据的基线，基于与均值和中位数的百分比差异、超出 Tukey 围栏上限以及使用 Z 分数来创建规则。我们发现，制定有效规则取决于仔细调整参数：均值和中位数差异的百分比、Tukey 围栏的乘数以及 Z 分数的阈值。为了确定哪个规则表现最好，我们使用了漏报率、假漏报率、假发现率和误报率。

在接下来的两章中，我们将介绍使用scikit-learn进行 Python 中的机器学习，并在第十一章《机器学习异常检测》中，我们将回顾这个场景，通过机器学习进行异常检测。

练习

完成以下练习以巩固本章所涵盖的概念：

1. 对 2018 年 12 月的数据进行模拟，并将结果保存到新的日志文件中，而无需重新创建用户基础。确保运行python3 simulate.py -h来查看命令行参数。将种子值设置为27。这些数据将用于剩余的练习。

2. 找出每个 IP 地址的唯一用户名数、尝试次数、成功次数、失败次数，以及成功/失败率，使用的是从练习1中模拟的数据。

3. 创建两个子图，左侧是失败次数与尝试次数的关系，右侧是失败率与唯一用户名数的关系。为结果图绘制决策边界。确保按是否为黑客 IP 地址为每个数据点着色。

4. 构建一个基于中位数百分比差异的规则，标记一个 IP 地址，如果其失败次数和尝试次数都为各自中位数的五倍，或者唯一用户名数为中位数的五倍。确保使用一个小时的时间窗口。记得使用get_baselines()函数计算基线所需的指标。

5. 使用本章中的evaluate()和classification_stats()函数计算指标，以评估这些规则的执行效果。

进一步阅读

查看以下资源，获取更多关于本章所讨论主题的信息：

• 引导法温和入门: machinelearningmastery.com/a-gentle-introduction-to-the-bootstrap-method/

• 引导法简介: towardsdatascience.com/an-introduction-to-the-bootstrap-method-58bcb51b4d60

• 向哈希算法中添加盐：存储密码的更好方式：auth0.com/blog/adding-salt-to-hashing-a-better-way-to-store-passwords/

• 暴力破解攻击：en.wikipedia.org/wiki/Brute-force_attack

• 分类准确率不足：更多可以使用的性能衡量标准：machinelearningmastery.com/classification-accuracy-is-not-enough-more-performance-measures-you-can-use/

• 字典攻击：en.wikipedia.org/wiki/Dictionary_attack

• 尼尔森规则：en.wikipedia.org/wiki/Nelson_rules

• 离线密码破解：攻击与最佳防御：www.alpinesecurity.com/blog/offline-password-cracking-the-attack-and-the-best-defense-against-it

• 泊松点过程：en.wikipedia.org/wiki/Poisson_point_process

• 精准度与召回率：en.wikipedia.org/wiki/Precision_and_recall

• Python 中的概率分布：www.datacamp.com/community/tutorials/probability-distributions-python

• 彩虹表：你密码的噩梦：www.lifewire.com/rainbow-tables-your-passwords-worst-nightmare-2487288

• RFC 1597（私有互联网的地址分配）：www.faqs.org/rfcs/rfc1597.html

• 抽样技术：towardsdatascience.com/sampling-techniques-a4e34111d808

• 修剪估计量：en.wikipedia.org/wiki/Trimmed_estimator

• 西电规则：en.wikipedia.org/wiki/Western_Electric_rules