您现在的位置是：首页 >技术交流 >《量化绿皮书》Chapter 7 Algorithms and Numerical Methods 算法与数值方法网站首页技术交流

《A Practical Guide To Quantitative Finance Interviews》，被称为量化绿皮书，是经典的量化求职刷题书籍之一，包含以下七章：

Chapter 1 General Principles 通用技巧

Chapter 2 Brain Teasers 脑筋急转弯

Chapter 3 Calculus and Linear Algebra 微积分与线性代数

Chapter 4 Probability Theory 概率论

Chapter 5 Stochastic Process and Stochastic Calculus 随机过程与随机微积分

Chapter 6 Finance 金融

Chapter 7 Algorithms and Numerical Methods 算法与数值方法

---

目录

 

7.1 Algorithms 算法

在编程中，算法复杂性的分析通常使用渐近分析，研究当输入的数量$n\to \infty$时的运行时间$T(n)$，即基本操作（如加法、乘法和比较）次数。

算法复杂度中的三种最重要的符号：

• $O$符号（描述渐近上界 asymptotic upper bound）： O ( g ( n ) ) = { f ( n ) : 存在正的常数 c , n 0 , 使得 0 ≤ f ( n ) ≤ c g ( n ) , ∀ n ≥ n 0 } O(g(n))={f(n):存在正的常数c,n_0,使得0le f(n)le cg(n),forall nge n_0} O(g(n))={f(n):存在正的常数c,n0​,使得0≤f(n)≤cg(n),∀n≥n0​}
• $\Omega$符号（描述渐近下界 asymptotic lower bound）： Ω ( g ( n ) ) = { f ( n ) : 存在正的常数 c , n 0 , 使得 0 ≤ c g ( n ) ≤ f ( n ) , ∀ n ≥ n 0 } Omega(g(n))={f(n):存在正的常数c,n_0,使得0le cg(n)le f(n),forall nge n_0} Ω(g(n))={f(n):存在正的常数c,n0​,使得0≤cg(n)≤f(n),∀n≥n0​}
• $\Theta$符号（描述渐近紧界 asymptotic tight bound）： Θ ( g ( n ) ) = { f ( n ) : 存在正的常数 c 1 , c 2 , n 0 , 使得 c 1 g ( n ) ≤ f ( n ) ≤ c 2 g ( n ) , ∀ n ≥ n 0 } Theta(g(n))={f(n):存在正的常数c_1,c_2,n_0,使得c_1g(n)le f(n)le c_2g(n),forall nge n_0} Θ(g(n))={f(n):存在正的常数c1​,c2​,n0​,使得c1​g(n)≤f(n)≤c2​g(n),∀n≥n0​}

算法复杂度中的两个概念：

• 最坏情况运行时间 worst-case running time：$W(n)$
  • 任意$n$个输入运行时间的上界
• 平均情况运行时间 average-case running time：$A(n)$
  • 随机选择$n$个输入的期望运行时间

对于许多算法，$W(n)$和$A(n)$具有相同的$O(g(n))$。但是，在某些问题中它们很可能是不同的，它们的相对重要性往往取决于具体问题。

**分治 divide-and-conquer：**一个有$n$个输入的问题通常可以分成$a$个子问题，每个子问题有$n/b$个输入。

• 如果需要$f(n)$次基本操作才能将问题划分为子问题并合并子问题的解，则运行时间可以表示为递归式$T(n)=aT(n/b)+f(n),a\ge 1,b\ge 1,f(n)\ge 0$

**主定理 master theorem：**求解递归方程$T(n)=aT(n/b)+f(n)$的紧界。

• 如果$f(n)=O(n^{{\log }_{b}a-ϵ})$，对某个常数$ϵ>0$，说明$f(n)$增长速度慢于$n^{{\log }_{b}a}$，因此$T(n)=\Theta (n^{{\log }_{b}a})$
• 如果$f(n)=\Theta (n^{{\log }_{b}a}{\log }^{k}n)$，对$k\ge 0$，说明$f(n)$和$n^{{\log }_{b}a}$与增长速度相似，因此$T(n)=\Theta (n^{{\log }_{b}a}{\log }^{k+1}n)$
• 如果$f(n)=\Omega (n^{{\log }_{b}a+ϵ})$，对某个常数$ϵ>0$，且$af(n/b)\le cf(n)$，对某个常数$c<1$，说明$f(n)$增长速度快于$n^{{\log }_{b}a}$，因此$T(n)=\Theta (f(n))$

主定理的应用：

• 二分查找 binary search：$T(n)=T(\frac{n}{2})+\Theta (1)$，$a=1,b=2,f(n)=\Theta (1),\Theta (n^{{\log }_{2}1})=\Theta (1)$，符合情况2，$T(n)=\Theta (\log n)$
• 归并排序 merge sort：$T(n)=2T(\frac{n}{2})+\Theta (n)$，$a=2,b=2,f(n)=\Theta (n),\Theta (n^{{\log }_{2}2})=\Theta (n)$，符合情况2，$T(n)=\Theta (n\log n)$
• 二叉树遍历 binary tree traversal：$T(n)=2T(\frac{n}{2})+O(1)$，$a=2,b=2,f(n)=O(1),\Theta (n^{{\log }_{2}2})=\Theta (n)$，符合情况1，$T(n)=\Theta (n)$

7.1.1 Number swap 交换数字

如何在不使用额外存储空间的情况下交换两个整数$i$和$j$?

How do you swap two integers, $i$ and $j$, without using additional storage space?

答案：

• 储存$i$和$j$之和，再减去$j$以提取$i$的值并赋给$j$，最后减去$j$（此时已是$i$的值）以提取$j$的值并赋给$i$

void swap(int &i, int &j){
    i = i + j;	// store the sum of i and j
    j = i - j;	// change j to i's value
    i = i - j;	// change i to j's value
}

• 按位异或操作 (XOR)，利用$x$^$x=0$和$0$$x=x$

void swap(int &i, int &j){
    i = i ^ j;
    j = i ^ j;	// j = (i ^ j) ^ j = i
    i = j ^ i;	// i = i ^ (i ^ j) = j
}

7.1.2 Unique elements 唯一元素

如果给你一个排序数组，你能写一些代码从数组中提取唯一的元素吗?例如，如果数组为[1,1,3,3,5,5,5,9,9,9,9]，则唯一元素应为[1,3,5,9]。

If you are given a sorted array, can you write some code to extract the unique elements from the array? For example, if the array is [1, 1, 3, 3, 3, 5, 5, 5, 9, 9, 9, 9], the unique elements should be [1, 3, 5, 9].

答案：

将数组按顺序排列，与前一元素不同的即为唯一元素

template <class T> vector<T> unique(T a[], int n){
    vector<T> vec;	// vector used to avoid resizing problem
    vec.reserve(n);	// reserver to avoid reallocation
    vec.push_back(a[0]);
    for(int i=1; i<n; ++1){
        if(a[i] != a[i-1])
            vec.push_back(a[i]);
    }
    return vec;
}

7.1.3 Horner’s algorithm 霍纳算法

写一个算法来计算$y=A_0+A_{1}x+A_2{x}^2+A_3{x}^3+\dots +A_n{x}^n$。

Write an algorithm to compute $y=A_0+A_{1}x+A_2{x}^2+A_3{x}^3+\dots +A_n{x}^n$.

答案：

计算每项并相加需要$O(n^2)$次乘法，霍纳算法可以将乘法次数减少到$O(n)$：

$y=((((A_{n}x+A_{n-1})x+A_{n-2})x+\dots +A_2)x+A_1)x+A_0$

逐项计算$B_n=A_n,B_{n-1}=B_{n}x+A_{n-1},\dots ,B_0=B_{1}x+A_0=y$

7.1.4 Moving average 移动平均

给定长度为m的大型数组a，您是否可以开发一种有效的算法来构建包含原始数组的n元素移动平均值的另一个数组（$B_1,\dots ,B_{n-1}=NA,B_i=(A_{i-n+1}+A_{i-n+2}+\dots +A_i)/n,\forall i=n,\dots ,m$）？

Given a large array A of length m, can you develop an efficient algorithm to build another array containing the n-element moving average of the original array ($B_1,\dots ,B_{n-1}=NA,B_i=(A_{i-n+1}+A_{i-n+2}+\dots +A_i)/n,\forall i=n,\dots ,m$)?

答案：

计算新的移动平均值时，可以复用前一值。前一移动平均值乘n，再减去前一移动平均值的第一项并加上新一项，再除n即可得到新的移动平均值。

S = A[1] + … + A[n]; B[n] = S/n;
for(i = n+1 to m){S = S - A[i-n] + A[i]; B[i] = S/n;}

7.1.5 Sorting algorithm 排序算法

你能解释一下对n个不同的值A1,…,An进行排序的三种排序算法，并分析每种算法的复杂度吗？

Could you explain three sorting algorithms to sort n distinct values A1,…,An and analyze the complexity of each algorithm?

答案：

• 插入排序（Insertion Sort）：将数组分为两个区间：左侧为有序区间，右侧为无序区间。每趟从无序区间取出一个元素，然后将其插入到有序区间的适当位置。

  • 时间复杂度：$O(n^2)$，需要$\frac{n(n-1)}{2}$次元素比较。
  • 空间复杂度：$O(1)$，原地排序算法，只用到指针变量$i$、$j$、无序区间的第1个元素$base$。

  

• 归并排序（Merge Sort）：采用经典的分治策略，先递归地将当前数组平均分成两半，然后将有序数组两两合并，最终合并成一个有序数组。

  • 时间复杂度：$O(n\ast logn)$，外层循环$\log n$量级，子算法 merge 为$n$量级。
  • 空间复杂度：$O(n)$，归并排序方法需要用到与参加排序的数组同样大小的辅助空间。

  

• 快速排序（Quick Sort）：采用经典的分治策略，选择数组中某个元素作为基准数，通过一趟排序将数组分为独立的两个子数组，一个子数组中所有元素值都比基准数小，另一个子数组中所有元素值都比基准数大。然后再按照同样的方式递归的对两个子数组分别进行快速排序，以达到整个数组有序。

  • 时间复杂度：$O(n\ast \log n)$。

  • 空间复杂度：$O(n)$，无论快速排序算法递归与否，排序过程中都需要用到堆栈或其他结构的辅助空间来存放当前待排序数组的首、尾位置。最坏的情况下，空间复杂度为$O(n)$。

  • 哨兵划分：选取一个基准数，将数组中比基准数大的元素移动到基准数右侧，比他小的元素移动到基准数左侧。

    1. 从数组中找到一个基准数$pivot$（以数组第1个元素作为基准数时$pivot=nums$）

    2. 使用指针 iii 指向数组开始位置，指针 jjj 指向数组末尾位置。

    3. 从数组中找到一个基准数$pivot$（默认为第一个元素$nums[low]$），交换到数组开始的位置。

    4. 指针$i$指向数组开始的位置，指针$j$指向数组末尾的位置。

       • 从右向左移动指针$j$，找到第一个小于基准数的元素位置。
       • 从左向右移动指针$i$，找到第一个大于基准数的元素位置。
       • 交换指针$i$和指针$j$所指位置的两个元素。

    5. 重复上述的三个步骤，直到指针$i$和指针$j$相遇时停止，最后将基准数与指针$i$所指位置的元素交换，将基准数放到两个子数组交界的位置。

  • 递归分解：完成哨兵划分之后，对划分好的左右子数组分别进行递归排序。

    1. 按照基准数的位置将数组拆分为左右两个子数组。
       1. 对每个子数组分别重复「哨兵划分」和「递归分解」，直到各个子数组只有1个元素，排序结束。

7.1.6 Random permutation 随机排列

A. 如果你有一个随机数生成器，可以从离散或连续均匀分布中生成随机数，你如何洗牌一副52张牌，使每种排列都是等可能的?

A. If you have a random number generator that can generate random numbers from either discrete or continuous uniform distributions, how do you shuffle a deck of 52 cards so that every permutation is equally likely?

答案：

给每张卡片分配一个随机数，然后按照分配的随机数对卡片进行排序（如果使用连续均匀分布，理论上任意两个随机数相等的概率为0）.通过对称，在$n!$个可能的顺序都是等可能的。复杂度由排序步骤决定，所以运行时间为$\Theta (n\log n)$

对于较大的n，我们可能需要使用一种更快的算法，称为Knuth shuffle。对于n个元素$A[1],\dots ,A[n]$，使用以下循环生成随机排列：

for(i=1 to n)swap(A[i], A[Random(i,n)])	//其中Random(i,n)是i和n之间的离散均匀分布中的随机数

复杂度为$\Theta (n)$，并且具有直观的解释。

• 在第一步中，n张牌中的每一张都有相同的概率被选为第一张牌，因为卡号是从1到n之间的离散均匀分布中选择的。
• 在第二步中，剩下的n -1张牌元素中的每一个都有相同的概率被选为第二张牌
• ……

所以每个有序序列都有1/n!的概率。

B. 你有一个由字符组成的文件。文件中的字符可以顺序读取，但文件的长度是未知的。如何选择一个字符，使文件中的每个字符被选中的概率相等?

B. You have a file consisting of characters. The characters in the file can be read sequentially, but the length of the file is unknown. How do you pick a character so that every character in the file has equal probability of being chosen?

答案：

从选择第一个字符开始：

• 如果有第二个字符，我们保留第一个字符的概率为1/2，并用概率为1/2的第二个字符替换选择。

• 如果有第三个字符，我们**保留选择（保留第一个字符或替换为第二个字符）**的概率为2/3，并用概率为1/3的第三个字符替换选择。

• 同样的过程一直持续到最后一个字符。

即设$C_n$是我们扫描n个字符后选择的字符，并且存在第$n+1$个字符，保留选择的概率为$\frac{n}{n+1}$，替换为第$n+1$个字符的概率为$\frac{1}{n+1}$。

使用简单的归纳法，我们可以很容易地证明，如果有m个字符，每个字符被选中的概率是1/m。

7.1.7 Search algorithm 搜索算法

A. 开发一种算法，使用不超过3n/2的比较来找到n个数的最小值和最大值。

A. Develop an algorithm to find both the minimum and the maximum of n numbers using no more than 3n/2 comparisons.

答案：

对于包含n个数字的未排序数组，需要进行最多n -1次比较才能确定数组的最小值或最大值。但可以通过最多3n/2次比较同时确定最小值和最大值。

• 如果我们将元素分成n/2对，比较每对中的元素，将较小的元素放入组A中，将较大的元素放入组B中，需要n/2次比较。
• 此时整个数组的最小值一定在组A中而最大值一定在组B中，因此只需要找到A中的最小值和B中的最大值，分别需要进行n/2-1次比较。

所以总比较次数最多是3n/2次。

B. 给你一组数字。从数组的开始到某个位置，所有元素都为零；在该位置之后，所有元素都是非零的。如果不知道数组的大小，如何找到第一个非零元素的位置?

B. You are given an array of numbers. From the beginning of the array to some position, all elements are zero; after that position, all elements are nonzero. If you don’t know the size of the array, how do you find the position of the first nonzero element?

答案：

从第一个元素开始：

• 如果它是零，我们检查第二个元素
• 如果第二个元素是零，我们检查第四个元素
• ……
• 重复这个过程，直到第$2^i$个元素非零

从第$2^i$个元素往回走：

• 检查第$\frac{2^i+2^{i-1}}{2}$个元素，如果它是零，搜索范围限制在第$\frac{2^i+2^{i-1}}{2}$个元素和第$2^i$个元素之间；否则搜索范围限制在第$2^{i-1}$个元素和第$\frac{2^i+2^{i-1}}{2}$个元素之间
• ……
• 重复这个过程，每次搜索范围减半

这种方法基本上是一个二分搜索。如果第一个非零元素位于位置n，则算法复杂度为$\Theta (\log n)$。

C. 你有一个正方形的数字网格。每一行的数字从左到右递增。每列中的数字从上到下递增。设计一个算法从网格中找到一个给定的数字。你的算法的复杂度是多少？

C. You have a square grid of numbers. The numbers in each row increase from left to right. The numbers in each column increase from top to bottom. Design an algorithm to find a given number from the grid. What is the complexity of your algorithm?

答案：

设数字网络是一个$n\times n$的矩阵$A$，$x$是我们要在网格中找到的数字。

• 在最后一列从上到下开始搜索：$A_{1,n},\dots ,A_{n,n}$

  • 如果找到数字，则停止搜索
  • 如果$A_{n,n}<x$，则$x$不在网格中，也停止搜索
  • 如果$A_{i,n}<x<A_{i+1,n}$，则排除第 $1,\dots ,i$行

• 在第$i+1$行从右到左开始搜索：$A_{i+1,n},\dots ,A_{i+1,1}$

  • 如果找到数字，则停止搜索
  • 如果$A_{i+1,1}>x$，则$x$不在网格中，也停止搜索
  • 如果$A_{i+1,j}<x<A_{i+1,j+1}$，则排除第$j+1,\dots ,n$列

• 在第$j$列从上到下沿着$A_{i+1,j}$到$A_{n,j}$搜索

  • 找到数字或$x$不在网格中，则停止搜索
  • 如果$A_{k,j}<x<A_{k+1,j}$，则排除第 $i+1,\dots ,k$行

• 在第$k+1$行从右到左开始搜索

• ……

使用这种算法，搜索最多需要2n步（从上到下的等价于一条线n，从右到左的等价于一条线n），复杂度为$O(n)$

7.1.8 Fibonacci numbers 斐波那契数

考虑以下生成斐波那契数的c++程序：

int Fibonacci(int n)
{
    if (n <= 0)
        return 0;
    else if (n == 1)
        return 1;
    else
        return Fibonacci(n-1) + Fibonacci(n-2)；
}

如果对于某个较大的n，计算Fibonacci(n)需要100秒，那么计算Fibonacci(n+1)需要多长时间，精确到秒？这个算法有效吗？如何计算斐波那契数列?

Consider the following C++ program for producing Fibonacci numbers:

int Fibonacci(int n)
{
    if (n <= 0)
        return 0;
    else if (n == 1)
        return 1;
    else
        return Fibonacci(n-1) + Fibonacci(n-2)；
}

If for some large n, it takes 100 seconds to compute Fibonacci(n), how long will it take to compute Fibonacci(n+1), to the nearest second? Is this algorithm efficient? How would you calculate Fibonacci numbers?

答案：

斐波那契数列定义为以下递归式：$F_0=0,F_1=1,F_n=F_{n-1}+F_{n-2},\forall n\ge 2$

通式：$F_n=\frac{(1+\sqrt{5}{)}^n-(1-\sqrt{5}{)}^n}{2^n\sqrt{5}}$

从c++程序可以推出运行时间序列：$T(0)=1,T(1)=1,T(n)=T(n-1)+T(n-2)+1$，因此运行时间也和斐波那契数列成正比，通式也近似。

对于某个较大的n，$(1-\sqrt{5}{)}^n\to 0$，因此$\frac{T(n+1)}{T(n)}\approx \frac{(1+\sqrt{5})}{2}$，$T(n+1)\approx \frac{(1+\sqrt{5})}{2}T(n)\approx 162$秒。

该递归算法的指数复杂度为$\Theta ((\frac{\sqrt{5}+1}{2}{)}^n)$，效率很低。原因是它不能有效地利用斐波那契数列中较小的斐波那契数的信息。如果按定义顺序计算$F_0,F_1,\dots ,F_n$，则运行时间的复杂度为$\Theta (n)$。

递归平方算法可以进一步将复杂度降低到$\Theta (\log n)$。

$\left[\begin{array}{cc}{F}_{n+1}& F_n\\ F_n& F_{n-1}\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\times \left[\begin{array}{cc}{F}_n& F_{n-1}\\ F_{n-1}& F_{n-2}\end{array}\right],\left[\begin{array}{cc}{F}_2& F_1\\ F_1& F_0\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$

归纳法可得：$\left[\begin{array}{cc}{F}_{n+1}& F_n\\ F_n& F_{n-1}\end{array}\right]={\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^n$

令$A=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$，使用分治法计算得$A^n=\{\begin{array}{ll}{A}^{n/2}\times A^{n/2},& 如果n为偶数\\ A^{(n-1)/2}\times A^{(n-1)/2}\times A,& 如果n为奇数\end{array}$

两个2x2矩阵的乘法复杂度为$\Theta (1)$，则$T(n)=T(n/2)+\Theta (1)$。应用主定理得到$T(n)=\Theta (\log n)$

7.1.9 Maximum contiguous subarray 最大连续子数组

假设你有一个长度为n的一维数组A，其中包含正数和负数。设计一种算法，求任意连续子数组$A[i,j]$的最大和$V(i,j)={\sum }_{x=i}^{j}A[x],1\le i\le j\le n$

Suppose you have a one-dimensional array A with length n that contains both positive and negative numbers. Design an algorithm to find the maximum sum of any contiguous j subarray $A[i,j]ofA:V(i,j)={\sum }_{x=i}^{j}A[x],1\le i\le j\le n$.

答案：

• 顺序计算，$O(n^2)$

$V(i,i)=A[i]$ when $j=i$ and $V(i,j)={\sum }_{x=i}^{j}A[x]=V(i,j-1)+A[j]$ when $j>i$

• 分治法，$O(n)$

$T(i)={\sum }_{x=1}^{i}A[x],T(0)=0$，则$V(i,j)=T(j)-T(i-1),\forall 1\le i\le j\le n$

显然对于任意固定$j$，当$T(i-1)$最小时，$V(i,j)$最大，所以以$j$结尾的最大子数组是$V_{\max }=T(j)-T_{\min },T_{\min }=\min (T(1),\dots ,T(j-1))$

在$j$增加时跟踪并更新$V_{\max }$和$T_{\min }$，可以开发以下算法：

T = A[1]; V_max = A[1]; T_min = min(0,T)
For j = 2 to n
{
	T = T + a[j];
	If T - T_min > V_max, then V_max = T - T_min;
	If T < T_min, then T_min = T;
}
Return V_max;

相应的c++函数，返回$V_{\max }$并给出数组及其长度的索引$i$和$j$：

double maxSubarray(double A[], int len, int &i, int &j)
{
    double T=A[0], Vmax=A[0];
    double Tmin = min(0.0, T);
    for(int k=1; k<len; Vmax)
    {
        T+=A[k];
        if (T-Tmin > Vmax) {Vmax = T-T_min; j=k;}
        if (T < Tmin) {Tmin = T; i = (k+1<j)? (k+1):j;}
    }
    return Vmax;
}


double A[] = {1.0, 2.0, -5.0, 4.0, -3.0, 2.0, 6.0, -5.0, -1.0};
int i = 0, j = 0;
double Vmax = maxSubarray(A, sizeof(a)/sizeof(A[1]), i, j);
// Vmax = 9, i = 3, j = 6, [4.0, -3.0, 2.0, 6.0]

7.2 The Power of Two 2的幂

7.2.1 Power of 2 2的幂

如何判断一个整数是否是2的幂？

How do you determine whether an integer is a power of 2?

答案：

任何整数$x=2^n(n\ge 0)$都有一位为1（从右数第$n+1$位），$8=2^3$表示为$0\dots 01000$。

同样，$2^n-1$从右数$n$位都为1，$7=2^3-1$表示为$0\dots 01000$。

所以$2^n$和$2^n-1$不共享任何位。因此，$x\&(x-1)=0$，其中$\&$是位与运算符，是识别整数x是否是2的幂的简单方法。

7.2.2 Multiplication by 7 乘7

给出一个快速的方法乘以一个整数7不使用乘法(*)操作符？

Give a fast way to multiply an integer by 7 without using the multiplication () operator?*

答案：

$(x<<3)-x$，其中$<<$是左移操作符，$x<<3$等价于$x\ast 8$

7.2.3 Probability simulation 概率模拟

给你一枚均匀的硬币。你能否设计一个使用公平硬币的简单游戏，使你获胜的概率为p， 0 < p < l?

You are given a fair coin. Can you design a simple game using the fair coin so that your probability of winning is p, 0 < p < 1?

答案：

将$p\in (0,1)$表示为二进制数，二进制数的每一位都可以用一个公平硬币来模拟。

$p=0.p_1{p}_2\dots p_n=p_1{2}^{-1}+p_2{2}^{-2}+\dots +p_n{2}^{-n},p_i\in (0,1),\forall i=1,2,\dots ,n$

抛均匀硬币，把正面数为1，反面数为0。设 s i ∈ { 0 , 1 } s_iin{0,1} si​∈{0,1}为从第$i$次投掷的结果，每次抛硬币后比较$p_i$和$s_i$。

• 如果$s_i<p_i$，我获胜，抛硬币停止
• 如果$s_i>p_i$，我输，抛硬币停止
• 如果$s_i=p_i$，继续抛硬币

例如$p=1/4=0\ast 2^{-1}+1\ast 2^{-2}=0.01$，只有抛出00才会获胜，01，10，11均输，获胜概率为1/4。

7.2.4 Poisonous wine 毒酒

你为生日聚会准备了1000瓶酒。派对开始前20小时，酒庄给你发了一条紧急消息，说有一瓶酒中毒了。你刚好有10只实验老鼠可以用来测试一瓶酒是否有毒。这种毒药非常强，任何剂量都能在18小时内杀死一只老鼠。但在死亡前18小时没有其他症状。在派对开始前，你能用这10只老鼠找到有毒的瓶子吗？

You’ve got 1000 bottles of wines for a birthday party. Twenty hours before the party, the winery sent you an urgent message that one bottle of wine was poisoned. You happen to have 10 lab mice that can be used to test whether a bottle of wine is poisonous.The poison is so strong that any amount will kill a mouse in exactly 18 hours. But before the death on the 18th hour, there are no other symptoms. Is there a sure way that you can find the poisoned bottle using the 10 mice before the party?

答案：

所有1 ~ 1000之间的整数都可以用10位二进制表示。例如，1000可以标记为1111101000。

• 让老鼠1从第1位（最低位）为1的瓶子中啜一口；
• 让老鼠2从第2位为1的瓶子里啜一口；
• …
• 让老鼠10从第10位（最高位）为1的瓶子中啜一口。

18小时后，如果我们把老鼠从最高位到最低位排成一排，把活老鼠当作0，把死老鼠当作1，我们就可以很容易地追溯毒瓶的标签。

例如，如果第六只、第七只和第九只老鼠都死了，而其他老鼠都还活着，则排列顺序为0101100000，有毒瓶子的标签为$2^8+2^6+2^5=352$。

n只老鼠可以判断$2^n$个瓶子。

7.3 Numerical Methods 数值方法

7.3.1 Monte Carlo simulation 蒙特卡洛模拟

蒙特卡罗模拟是一种使用具有适当概率的随机数作为输入，对确定性模型进行迭代评估的方法。对于衍生品定价，它以对应于标的随机过程的概率模拟大量标的资产的价格路径（通常在风险中性度量下），计算每条路径下的衍生品的折现回报，并将折现回报平均以得到衍生品价格。蒙特卡罗模拟的有效性依赖于大数定律。

如果收益只依赖于标的资产的最终价值，蒙特卡罗模拟可以用来估计衍生品的价格，如果收益也依赖于路径，蒙特卡罗模拟也可以用来估计价格。然而，它不能直接应用于美式期权或任何其他具有早期行权期权的衍生品。

A. 解释你如何使用蒙特卡罗模拟来为欧式看涨期权定价？

A. Explain how you can use Monte Carlo simulation to price a European call option?

答案：

如果我们假设股票价格遵循几何布朗运动，我们可以模拟可能的股票价格路径。

把$t$和$T$之间的时间分成$N$个等间隔的时间步长，即$\Delta t=\frac{T-t}{N},t_i=t+\Delta t\times i$，然后用方程$S_i=S_{i-1}{e}^{(r-{\sigma }^2/2)(\Delta t)+\sigma \sqrt{\Delta t}{ϵ}_i}$模拟风险中性概率下的股价路径（${ϵ}_i$为标准正态分布中的IID随机变量）。假设我们模拟$M$条路径，每个路径产生一个在到期日$T$的股票价格$S_{T,k}$。

欧式看涨期权的估计价格是预期收益的现值：$C=e^{-r(T-t)}\frac{{\sum }_{k=1}^M\max (S_{T,k}-K,0)}{M}$

B. 如果您的计算机只能生成遵循0到1之间连续均匀分布的随机变量，您如何生成遵循$N(\mu ,{\sigma }^2)$的随机变量？

B. How do you generate random variables that follow $N(\mu ,{\sigma }^2)$ if your computer can only generate random variables that follow continuous uniform distribution between 0 and 1?

答案：

这个问题的解决方法可以分为两个步骤：

1. 利用逆转换法（inverse transform method）和接受-拒绝法（acceptance-rejection method）从均匀随机数生成器中生成$x\sim N(0,1)$的随机变量。
2. 将$x$缩放到$\mu +\sigma x$，生成遵循$N(\mu ,{\sigma }^2)$的随机变量

逆转换法（inverse transform method）：对于任意具有累积密度函数$F(U=F(X))$的连续随机变量$X$，随机变量$X$可以定义为$U$的反函数$X=F^{-1}(U),0\le U\le 1$，很明显是一个一对一函数。因此，任何连续随机变量都可以通过以下过程生成：

• 从标准均匀分布中生成一个随机数u。
• 计算值$x$，使$u=F(x)$作为F描述的分布中的随机数（$F^{-1}(U)$必须是可计算的）。

对于标准正态分布，$U=F(X)={\int }_{-\infty }^X\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}dx$，逆函数没有解析解。我们用欧拉法等数值积分法，得到$X$到$U$的一对一映射作为常微分方程$F^{\prime }(x)=f(x)=\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}$的数值解。

接受-拒绝法（acceptance-rejection method）：一些随机变量有pdf $f(x)$但$F^{-1}(U)$没有解析解。在这些情况下，我们可以使用pdf $g(y)$和$Y=G^{-1}(U)$的随机变量$y$来帮助生成pdf $f(x)$的随机变量。

假设$M$是一个满足$\frac{f(y)}{g(y)}\le M,\forall y$的常数，

• 抽样：从$g(y)$生成随机变量$y$，从标准均匀分布[0,1]生成随机变量$v$。
• 接受/拒绝：如果$v\le \frac{f(y)}{Mg(y)}$，接受$x=y$；否则，重复抽样。

指数随机变量$g(x)=\lambda e^{-\lambda x},\lambda =1$有cdf $u=G(x)=1-e^{-x}$，逆函数有解析解$x=-\log (1-u)$，因此指数分布的随机变量可以方便地模拟。

对于标准正态分布，$f(X)=\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}$，$\frac{f(x)}{g(x)}=\sqrt{\frac{2}{\pi }}{e}^{x-x^2/2}<\sqrt{\frac{2}{\pi }}{e}^{-(x-1{)}^2/2+1/2}\le \sqrt{\frac{2}{\pi }}{e}^{1/2}\approx 1.32,\forall 0<x<\infty$，因此可以选取$M=1.32$并使用接受-拒绝法生成$x\sim N(0,1)$的随机变量，并将其缩放为$N(\mu ,{\sigma }^2)$的随机变量。

C. 你能解释一些方差减少技术来提高蒙特卡罗模拟的效率吗？

C. Can you explain a few variance reduction techniques to improve the efficiency of Monte Carlo simulation?

答案：

蒙特卡罗模拟的基本形式是IID随机变量$Y_1,Y_2,\dots ,Y_M$的均值$\overline{Y}=\frac{1}{M}{\sum }_{i=1}^M{Y}_i$。因为每个$Y_i$的期望值是无偏的，所以估计量$\overline{Y}$也是无偏的；如果$Var(Y)=\sigma$，则$Var(\overline{Y})=\sigma /\sqrt{M}$。如果$\sigma$很大，蒙特卡罗模拟的计算量很大，通常需要数千甚至数百万次模拟才能获得所需的精度。

根据具体问题的不同，已经应用了各种方法来减小方差：

对偶变量 Antithetic variable：对于每一系列${ϵ}_i$，计算其对应的$Y({ϵ}_1,\dots ,{ϵ}_n)$。然后把所有${ϵ}_i$的符号倒过来，并计算相应的$Y(-{ϵ}_1,\dots ,-{ϵ}_n)$。当$Y({ϵ}_1,\dots ,{ϵ}_n)$和$Y(-{ϵ}_1,\dots ,-{ϵ}_n)$负相关时，方差减小。

矩匹配 Moment matching：随机变量的特定样本可能不能很好地匹配总体分布。我们可以先绘制一大组样本，然后重新缩放样本，使样本的矩（最常用的是均值和方差）与期望的总体矩匹配。

控制变量 Control variate：如果我们想为导数$X$定价，并且有一个密切相关的导数$Y$有解析解，我们可以生成一系列随机数，并使用相同的随机序列为$X$和$Y$定价，从而得到$\hat{X}$和$\hat{Y}$，然后$X$可以估计为$\hat{X}+(Y-\hat{Y})$。本质上使用$(Y-\hat{Y})$来修正$X$的估计误差。

重要性抽样 Importance sampling：为了从分布$f(x)$中估计$h(x)$的期望值，我们可以从分布$g(x)$中绘制$x$，而不是从分布$f(x)$中绘制$x$。用蒙特卡罗模拟估计$\frac{h(x)f(x)}{g(x)}$的期望值：

$E_{f(x)}[h(x)]=\int h(x)f(x)dx=\int \frac{h(x)f(x)}{g(x)}g(x)dx=E_{g(x)}[\frac{h(x)f(x)}{g(x)}]$

如果$\frac{h(x)f(x)}{g(x)}$的方差小于$h(x)$，那么重要性抽样可以产生更有效的估计器。

这种方法可以用深度虚值期权作为例子来更好地解释：如果我们直接使用风险中性的$f(S_T)$作为分布，大多数模拟路径将产生$h(S_T)=0$，因此估计方差将很大。如果我们引入一个具有更宽跨度的分布$g(S_T)$（对于$S_T$来说尾部更宽），那么更多的模拟路径将具有正的$h(S_T)$。缩放因子$\frac{f(x)}{g(x)}$将保持估计量无偏，但具有更低的方差。

低差异序列 Low-discrepancy sequence：不使用随机样本，我们可以生成代表分布的“随机变量”的确定性序列。这种低差序列可以使收敛速度达到$1/M$。

D. 如果一个期权没有封闭形式的定价公式，你将如何估计它的delta和gamma？

D. If there is no closed-form pricing formula for an option, how would you estimate its delta and gamma?

答案：

如在A中，有或没有封闭形式定价公式的期权价格可以使用蒙特卡罗模拟推导出来，通过将当前的基础价格从$S$稍微改变为$S\pm \delta S$可以估计delta和gamma，其中$\delta S$是一个小的正值。

对所有三个起始价格$S-\delta S,S,S+\delta S$进行蒙特卡罗模拟，得到对应的期权价格$f(S-\delta S),f(S),f(S+\delta S)$。

$\Delta =\frac{\delta f}{\delta S}=\frac{f(S+\delta S)-f(S-\delta S)}{2\delta S}$

$\Gamma =\frac{(f(S+\delta S)-f(S))-(f(S)-f(S-\delta S))}{\delta S^2}$

为了减少方差，通常最好使用相同的随机数序列来估计$f(S-\delta S),f(S),f(S+\delta S)$。

E. 如何使用蒙特卡罗模拟来估计$\pi$？

E. How do you use Monte Carlo simulation to estimate $\pi$?

答案：

估计$\pi$的一种标准方法是在单位正方形中随机选择点（$x$和$y$是$0$到$1$之间的独立均匀随机变量），$\pi$为在1/4圆$x^2+y^2\le 1,x>0,y>0$内的点的比例的4倍。

7.3.2 Finite difference method 有限差分法

有限差分法是衍生品定价的另一种常用数值方法。它通过离散时间和标的证券的价格，数值求解微分方程来估计导数的价格。

将二阶非线性偏微分方程BSM方程转换为热扩散方程，这个新方程表示为$\tau$（到期时间）和$x$（标的证券价格的函数）的函数，是导数的一般微分方程。各种导数的区别在于边界条件。通过建立$x$和$\tau$的网格并使用边界条件，我们可以用有限差分方法递归地计算每个$x$和$\tau$处的$u$。

A. 你能简单解释一下有限差分法吗？

A. Can you briefly explain finite difference methods?

答案：

有限差分法在实际应用中分为显式差分法、隐式差分法和Crank-Nicolson法。

如图所示，将$\tau ,[0,T]$的范围划分为$N$个以增量$\Delta \tau =T/N$为单位的离散区间，将$x,[x_0,x_J]$的范围划分为$J$个以增量$\Delta x=(x_J-x_0)/J$为单位的离散区间，则时间$\tau$和空间$x$可以表示为${\tau }_n,n=1,\dots ，N$和$x_j,j=1,\dots ,J$的网格。

显式差分法：使用时间${\tau }_n$前向差分和$x_j$二阶中心差分

$\frac{\partial u}{\partial \tau }\approx \frac{u_j^{n+1}-u_j^n}{\Delta \tau }=\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{(\Delta x{)}^2}\approx \frac{{\partial }^{2}u}{\partial x^2}$

重新排列项，可以将$u_j^{n+1}$表示为$u_{j+1}^n,u_j^n,u_{j-1}^n$的线性组合$u_j^{n+1}=\alpha u_{j-1}^n+(1-2\alpha )u_j^n+\alpha u_{j+1}^n,其中\alpha =\Delta t/(\Delta x{)}^2$。此外，经常有边界条件$u_j^0,u_0^n,u_J^n,\forall n=1,\dots ,N;j=0,\dots ,J$。

结合边界条件和方程$u_j^{n+1}=\alpha u_{j-1}^n+(1-2\alpha )u_j^n+\alpha u_{j+1}^n,其中\alpha =\Delta t/(\Delta x{)}^2$，我们可以估计网格上所有的$u_j^n$。

隐式差分法：使用时间${\tau }_n$后向差分和$x_j$二阶中心差分

$\frac{\partial u}{\partial \tau }\approx \frac{u_j^{n+1}-u_j^n}{\Delta \tau }=\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{(\Delta x{)}^2}\approx \frac{{\partial }^{2}u}{\partial x^2}$

Crank-Nicolson法：使用时间$(t_n+t_{n+1})/2$的中心差分和$x_j$二阶中心差分

$\frac{\partial u}{\partial \tau }\approx \frac{u_j^{n+1}-u_j^n}{\Delta \tau }=\frac{1}{2}(\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{(\Delta x{)}^2}+\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{(\Delta x{)}^2})\approx \frac{{\partial }^{2}u}{\partial x^2}$

B. 如果用显式有限差分法求解抛物型偏微分方程，在时间维度上有太多步，还是在空间维度上有太多步，哪个更糟糕？

B. If you are solving a parabolic partial differential equation using the explicit finite difference method, is it worse to have too many steps in the time dimension or too many steps in the space dimension?

答案：

显式有限差分法中$u_j^{n+1}$的方程为$u_j^{n+1}=\alpha u_{j-1}^n+(1-2\alpha )u_j^n+\alpha u_{j+1}^n,其中\alpha =\Delta t/(\Delta x{)}^2$。

为了使显式有限差分法稳定，我们需要$1-2\alpha >0\Rightarrow \Delta t/(\Delta x{)}^2<1/2$。因此，较小的$\Delta y$（太多的时间步长）是理想的，但较小的$\Delta x$（太多的空间步长）可能使$\Delta t/(\Delta x{)}^2>1/2$，结果不稳定。

从这个意义上说，在空间维度上有太多步是更糟糕的。相比之下，隐式差分法总是稳定和收敛的。