作业1 用Monte Carlo方法进行概率和分位计算。

目的：提升统计软件的使用能力，熟练使用Monte Carlo方法进行概率计算，并能把Monte Carlo方法用于概率计算当成问题来进行科学研究，而不仅仅是完成计算，要当成一种方法进行研究。

基本要求：（1）最好选择一个需要计算较复杂统计量概率计算的问题，用Monte Carlo方法来完成；

（2）给出计算方法过程，输出一些计算结果；

（3）讨论Monte Carlo方法计算结果随着样本容量等量变化时计算误差的变化，并进行总结；

（4）与其他概率计算方法进行比较分析，给出Monte Carlo方法的优缺点。

（5）按发表论文格式整理。

**注：**除完成基本要求外，可以增加自己感兴趣的研究内容。

0 金融背景

在实际问题中，如果P_t 表示t时刻的价格，考虑以下不同收益率：

简单收益率 $R_t = \frac{P_t}{P_{t-1}} -1$

对数收益率 $r_t =ln(1+R_t)$

如果收益率服从正太分布， 一方面收益率不能取负值（当然在某些场景下，负收益率代表亏损也是有意义的）。另一方面，如果希望让分期越来越细，则随着分期次数越来越多，趋于无限，如果单纯将一期分为两期设置分期收益率为原来的一半，最后得到的总收益率会比原先更高（如下图所示）。

如果对数收益率服从正太分布，对于同一投资产品，对数收益率时序可加。

如果我们考察单一投资品在总共 T 期内的表现，那应该用对数收益率，而非算数收益率。**算术平均值不能正确的反应一个投资品的收益率。**比如一个投资品今年涨了 50%，明年跌了 50%，它的算数平均收益率为 0；但事实上，两年后该投资品亏损了最初资金的 25%。相反的，**对数收益率由于具备可加性，它的均值可以正确反映出该投资品的真实收益率。**比如这两年的对数收益率分别为 40.5% 和 -69.3%，平均值为 -28.77%，转换为百分比亏损就是 exp{-28.77%} - 1 = -25%。

对数收益率的时序可加性让我们能够使用另外两个利器：“中心极限定理”和“大数定律”。假设初始资金 $X_0$（假设等于 1），$ln(X_T) = ln(X_T/X_0)$ 就是整个 T 期的对数收益率。对数收益率的最大好处是它的可加性，把单期的对数收益率相加就得到整体的对数收益率。

如果能假设不同期是相互独立的，T 期对数收益率相加相当于 T 个独立的随机变量相加。由中心极限定理可知，它们的和逼近正态分布。由大数定律可知，$\frac{1}{T} ln\frac{X_T}{X_0}$ ，即单期对数收益率均值，随着 T 的增大一定会收敛于它的期望

我们以初始资金 $X_0$ 跑一个策略进行投资，最终是想让给定 T 期之后的 $X_T$ 越大越好，但我们不知道 $X_T$ 最终会收敛到什么值。但上面的分析说明只要 T 足够大，大数定律保证了 $X_T$ 的对数，即 $ln(X_T)$，一定会非常接近它的期望 $E[ln(X_T)]$，这就是用对数收益率的价值。

1 选取较复杂统计量（问题1）

假设对数收益率符合正太分布 $r \backsim N(\mu, \sigma^2)$，则收益率符合对数正太分布 $R \backsim Log-N(\mu, \sigma^2)$。收益率的分布的确也符合直觉，有少部分的投资者获取到单期高额收益率，而绝大多数投资者的收益率在0附近。

在金融学领域，一阶矩和二阶矩捕捉随机收益率的的长期（平均） 收益与风险; 三阶矩对极端风险的度量有重要指示 作用，四阶矩考察厚尾分布,与波动率的预测、估计 的有效性密切相关。

import math
import numpy as np
import matplotlib.pyplot as plt


u = 0   # 均值μ
sig = math.sqrt(1)  # 标准差δ
x = np.linspace(u - 3*sig, u + 3*sig, 50)   # 定义域
f = np.exp(-(x - u) ** 2 / (2 * sig ** 2)) / (math.sqrt(2*math.pi)*sig) # 定义曲线函数
y = np.exp(x) - 1
plt.figure(figsize=(20, 8))
plt.subplot(121)
plt.xlabel('X')
plt.ylabel('f')
plt.plot(x, f, "g", linewidth=2)    # 正太分布
plt.subplot(122)
plt.xlabel('Y')
plt.plot(y, f, "g", linewidth=2)    # 对数正太分布
plt.grid(True)  # 网格线
plt.show()  # 显示

现考虑**对数正太分布（收益率）**的偏度和峰度，设 $\overrightarrow T$：

$$\overrightarrow T= (\hat S, \hat K)$$

其中 $\hat S$ 表示变量的偏度，是统计数据分布偏斜方向和程度的度量，是统计数据分布非对称程度的数字特征。偏度(Skewness)亦称偏态、偏态系数。表征概率分布密度曲线相对于平均值不对称程度的特征数。直观看来就是密度函数曲线尾部的相对长度。

$$skewness = E[(\frac{X-\mu}{\sigma})^3] = \frac{1}{n}\sum_1^n(\frac{X_i-\mu}{\sigma})^3$$

对于偏度的样本估计（当$\mu$和$\sigma$ 不知道的情况下用 $\overline X$ 和 $S=\frac{\sum_1^n (X_i-\overline X)^2}{n-1}$ 代替）：

$$\hat {skewness} = \frac{\sum_1^n(X_i-\overline X)^3}{(n-1)S^3} = \frac{m_3}{m_2^{\frac{3}{2}}}$$

其中 $m_k$ 类似k阶中心距，但自由度少1，不妨记为类中心距，公式为：

$$m_k = \frac{\sum_1^n(X_i-\overline X)^k}{n-1}$$

在一般情形下，当统计数据为右偏分布时，S_k > 0，且 S_k 值越大，右偏程度越高；当统计数据为左偏分布时，S_k < 0，且 S_k 值越小，左偏程度越高。当统计数据为对称分布时，显然有 S_k = 0。

在实际计算中虽然可以采用定义式进行计算，但是需要多次遍历各个样本以计算均值和方差，当样本容量很大时耗时很大。所以常常利用1至3阶原点矩进行计算：

$\hat K$ 表示变量的峰度

$$kurtosis = E[(\frac{X-\mu}{\sigma})^4]$$

同样，在总体分布不清楚的情况下，考虑样本估计：

$$\hat{kurtosis} = \frac{\sum_1^n(X_i-\overline X)^4}{(n-1)S^4} = \frac{m_4}{m_2^2}$$

在更通常的情况下，峰度被定义为四阶累积量除以二阶累积量的平方，它等于四阶中心矩除以概率分布方差的平方再减去3，因此，本文也采用这种定义方式：

$$\hat{kurtosis} = \frac{m_4}{m_2^2} -3$$

2 模拟计算

若按照年化收益率近似5%的标准，取对数收益率符合正太分布 $N(0.05, 0.04)$，分布大致如图所示：左侧是对数收益率，右侧是简单收益率

从维基百科https://zh.wikipedia.org/wiki/对数正态分布上查询到该分布的偏度和峰度分别是：

偏度
峰度

带入 $\sigma^2= 0.04$ 得到

$$skewness = 0.6142947619866632 \\ kurtosis = 0.6783657771754363$$

将此值作为偏度和峰度的理论值。

但是希望计算偏度和峰度这两个统计量的固定区间上的概率及固定概率分位点却无法计算出来（过于复杂，而且没有现成的推导公式），因此后续的误差分析只能针对偏度和峰度来计算。

小样本模拟（N=1000）

对于给定t，计算$P\{\hat{skewness} \le t\}$，每次从$R_1,R_2...R_N$从取一个样本（这个样本维度为N=1000），并计算相应的统计量，如此反复，直到获取到samples = 1000个样本，得到如下表格：

对于结果进行统计描述，选取常见的统计量：均值，标准差，最大，最小，25%/50%/75%分位点值，分别为：

继而绘制偏度和峰度的分布：

从而计算偏度和峰度固定区间上的概率及固定概率分位点：为了简化问题，这里不妨取 $t_1 = t_2 = p_1 = p_2 = 0.5$，分位点与中位数相等，分别求出以下4个量：

$P\{skewness \le t_1\} = 0.146$

$F^{-1}_{skewness}(p) = 0.599008$

$P\{kurtosis\le t_2\} = 0.428$

$F^{-1}_{kurtosis}(p) = 0.560516$

计算误差：使用 skewness 和 kurtosis 的理论值，和实际得到的样本求出均方误差MSE及相对偏差：

$$MSE_{skewness} = \frac{1}{smaples}\sum_{i=1}^{samples}(skewness_i - skewness)^2 \\ \epsilon_{skewness} = |\frac{\overline{skewness}}{skewness}-1|$$

其中$\overline{skewness}$表示样本均值。kurtosis 计算方式同上.

P_skewness          0.146000
P_kurtosis          0.428000
F-1_skewness        0.599008
F-1_kurtosis        0.560516
MSE_skewness        0.011139
MSE_kurtosis        0.190482
epsilon_skewness    0.010996
epsilon_kurtosis    0.044816

扩大样本个数

将 N 与 samples同比例指数扩大，每次记录以上8个量，同时绘制两者的分布。此处取 N=samples = 1000, 2000，... 9000，分别得到结果如下图：

N	分布情况
1000
2000
3000
4000
5000
6000
7000
8000
9000

3 随着样本容量等量变化时计算误差的变化

绘制处以上统计结果随着样本容量变化的结果的统计图：

如图所示，随着样本容量逐渐提高，不论是MSE，还是epsilon都有显著的降低，这种降低对于kurtosis要远远明显于skewness，这可能因为，kurtosis是四阶的统计量，刻画的是波动率，对于数据量的需求更大。

4 Monte Carlo方法的优缺点

与其他概率计算方法进行比较分析，给出Monte Carlo方法的优缺点：

蒙特卡罗法的最大优点是:

1.方法的误差与问题的维数无关。

2.对于具有统计性质问题可以直接进行解决。

3.对于连续性的问题不必进行离散化处理

蒙特卡罗法的缺点则是:

1.对于确定性问题需要转化成随机性问题。

2.误差是概率误差。

3.通常需要较多的计算步数N.

蒙特卡罗法作为一种计算方法，是由美国数学家乌拉姆(Ulam , S. M.)与美籍匈牙利数学家冯·诺伊曼(von Neumann,J.)在20世纪40年代中叶，为研制核武器的需要而首先提出来的.实际上，该方法的基本思想早就被统计学家所采用了.例如，早在17世纪，人们就知道了依频数决定概率的方法。

附录

主体部分代码：

import numpy as np
import matplotlib.pyplot as plt
import math
import pandas as pd

mu = 0.05
sigma = 0.2
sigma2 = pow(sigma, 2)
e_sigma2 = math.exp(sigma2)
real = {
    'skewness': (e_sigma2 + 2) * math.sqrt(e_sigma2 - 1),
    'kurtosis': pow(e_sigma2, 4) + 2 * pow(e_sigma2, 3) + 3 * pow(e_sigma2, 2) - 6
}

static = pd.DataFrame(columns = [f'{b}_{a}' for b in ['P', 'F-1', 'MSE', 'epsilon'] for a in ['skewness', 'kurtosis']])

for N in range(1000, 10000, 1000):
    bins = int(N / 100)
    np.random.seed(0)
    samples = N

    df = pd.DataFrame(columns=[2, 3, 4, 'skewness', 'kurtosis'])
    for j in range(samples):
        # print(df.shape[0])
        r = np.random.normal(mu, sigma, N)
        R = np.exp(r) - 1
        R_mean = np.mean(R)
        # print(R - R_mean)

        m = {}
        for i in [2,3,4]:
            m[i] = np.sum(np.power(R - R_mean, i)) / (N-1)

        df.loc[df.shape[0]] = {**m, 
            "skewness": m[3] / math.sqrt(pow(m[2], 3)),
            "kurtosis": m[4] / pow(m[2], 2) - 3
        }

    df.rename(columns={a: f'm{a}' for a in df.columns if type(a) == int}, inplace=True)
    print(df)
    print(df.describe())
    df.describe().to_csv('df_describe.csv')
    df.to_csv('df.csv')

    threshold = 0.5
    dic = {}
    for key in ['skewness', 'kurtosis']:
        dic = {
            **dic, 
            f'P_{key}': sum(df[key] < threshold) / samples,
            f'F-1_{key}': df[key].quantile(threshold),
            f'MSE_{key}': np.mean(np.power(df[key] - real[key], 2)),
            f'epsilon_{key}': abs(df[key].mean() / real[key] - 1)
        }

    static.loc[static.shape[0]] = dic

    print(static.loc[0])

    plt.figure(figsize=(20, 8))
    plt.subplot(121)
    plt.xlabel('skewness')
    plt.ylabel('N')
    plt.hist(df['skewness'], bins)
    plt.subplot(122)
    plt.xlabel('kurtosis')
    plt.hist(df['kurtosis'], bins)
    plt.savefig(f'{N}.png')

static.to_csv('static.csv')