Python 中的 Monte Carlo 模拟使用使用 sobol 序列的准随机标准正态数给出了错误的值 [英] Monte Carlo simulations in Python using quasi random standard normal numbers using sobol sequences gives erroneous values

问题描述

我正在尝试使用准随机标准正态数执行蒙特卡罗模拟.我明白我们可以使用Sobol序列生成均匀数，然后使用概率积分变换将它们转换为标准正态数.我的代码给出了模拟资产路径的不切实际的值:

import sobol_seq将 numpy 导入为 np从 scipy.stats 导入规范def i4_sobol_generate_std_normal(dim_num, n, skip=1):"""生成多元标准正态准随机变量.参数:输入，整数dim_num，空间维度.输入，整数 n，要生成的点数.输入，整数 SKIP，要跳过的初始点数.输出，形状(n，dim_num)的真实np数组."""sobols = sobol_seq.i4_sobol_generate(dim_num, n, 跳过)法线 = norm.ppf(sobols)回归正常def GBM(Ttm, TradingDaysInAYear, NoOfPaths, UnderlyingPrice, RiskFreeRate, Volatility):dt = float(Ttm)/TradingDaysInAYear路径 = np.zeros((TradingDaysInAYear + 1, NoOfPaths), np.float64)路径[0] = 标的价格对于范围内的 t(1, TradingDaysInAYear + 1):rand = i4_sobol_generate_std_normal(1, NoOfPaths)lRand = []对于范围内的 i(len(rand)):a = rand[i][0]lRand.append(a)rand = np.array(lRand)路径[t] = 路径[t - 1] * np.exp((RiskFreeRate - 0.5 * 波动率** 2) * dt + 波动率* np.sqrt(dt) * rand)返回路径GBM(1, 252, 8, 100., 0.05, 0.5)数组([[1.00000000e+02, 1.00000000e+02, 1.00000000e+02, ...,1.00000000e+02, 1.00000000e+02, 1.00000000e+02],[9.99702425e+01, 1.02116774e+02, 9.78688323e+01, ...,1.00978615e+02, 9.64128959e+01, 9.72154915e+01],[9.99404939e+01, 1.04278354e+02, 9.57830834e+01, ...,1.01966807e+02, 9.29544649e+01, 9.45085180e+01],...,[9.28295879e+01, 1.88049044e+04, 4.58249200e-01, ...,1.14117599e+03, 1.08089096e-02, 8.58754653e-02],[9.28019642e+01, 1.92029616e+04, 4.48483141e-01, ...,1.15234371e+03, 1.04211828e-02, 8.34842557e-02],[9.27743486e+01, 1.96094448e+04, 4.38925214e-01, ...,1.16362072e+03, 1.00473641e-02, 8.11596295e-02]])

不应生成 8.11596295e-02 之类的值，因此我认为代码中有问题.如果我使用来自 numpy 库 rand = np.random.standard_normal(NoOfPaths) 的标准法线绘制，那么价格与 Black Scholes 的价格相匹配.因此我认为问题出在随机数生成器上.值 8.11596295e-02 指的是路径中的价格，价格不太可能从 100(初始价格)下降到 8.11596295e-02.

参考文献:

I'm trying to perform Monte Carlo Simulations using quasi-random standard normal numbers. I understand that we can use Sobol sequences to generate uniform numbers, and then use probability integral transform to convert them to standard normal numbers. My code gives unrealistic values of the simulated asset path:

import sobol_seq
import numpy as np
from scipy.stats import norm

def i4_sobol_generate_std_normal(dim_num, n, skip=1):
    """
    Generates multivariate standard normal quasi-random variables.
    Parameters:
      Input, integer dim_num, the spatial dimension.
      Input, integer n, the number of points to generate.
      Input, integer SKIP, the number of initial points to skip.
      Output, real np array of shape (n, dim_num).
    """

    sobols = sobol_seq.i4_sobol_generate(dim_num, n, skip)

    normals = norm.ppf(sobols)

    return normals

def GBM(Ttm, TradingDaysInAYear, NoOfPaths, UnderlyingPrice, RiskFreeRate, Volatility):
    dt = float(Ttm) / TradingDaysInAYear
    paths = np.zeros((TradingDaysInAYear + 1, NoOfPaths), np.float64)
    paths[0] = UnderlyingPrice
    for t in range(1, TradingDaysInAYear + 1):
        rand = i4_sobol_generate_std_normal(1, NoOfPaths)
        lRand = []
        for i in range(len(rand)):
            a = rand[i][0]
            lRand.append(a)
        rand = np.array(lRand)

        paths[t] = paths[t - 1] * np.exp((RiskFreeRate - 0.5 * Volatility ** 2) * dt + Volatility * np.sqrt(dt) * rand)
    return paths

GBM(1, 252, 8, 100., 0.05, 0.5)

array([[1.00000000e+02, 1.00000000e+02, 1.00000000e+02, ...,
        1.00000000e+02, 1.00000000e+02, 1.00000000e+02],
       [9.99702425e+01, 1.02116774e+02, 9.78688323e+01, ...,
        1.00978615e+02, 9.64128959e+01, 9.72154915e+01],
       [9.99404939e+01, 1.04278354e+02, 9.57830834e+01, ...,
        1.01966807e+02, 9.29544649e+01, 9.45085180e+01],
       ...,
       [9.28295879e+01, 1.88049044e+04, 4.58249200e-01, ...,
        1.14117599e+03, 1.08089096e-02, 8.58754653e-02],
       [9.28019642e+01, 1.92029616e+04, 4.48483141e-01, ...,
        1.15234371e+03, 1.04211828e-02, 8.34842557e-02],
       [9.27743486e+01, 1.96094448e+04, 4.38925214e-01, ...,
        1.16362072e+03, 1.00473641e-02, 8.11596295e-02]])

Values like 8.11596295e-02 should not be generated, hence I think there is something wrong in the code. If I use standard normal draws from the numpy library rand = np.random.standard_normal(NoOfPaths) then the price matches with the Black Scholes price. Hence I think the problem is with the random number generator. The value 8.11596295e-02 refers to a price in a path, and it's very unlikely that the price would come down from 100 (initial price) to 8.11596295e-02.

References: 1, 2, 3.

解决方案

It appears there is a bug in sobol_seq. Anaconda, python 3.7, 64bit, Windows 10 x64, installed sobol_seq via pip

pip install sobol_seq

# Name                    Version                   Build  Channel
sobol-seq                 0.1.2                    pypi_0    pypi

Simple code

print(sobol_seq.i4_sobol_generate(1, 1, 0))
print(sobol_seq.i4_sobol_generate(1, 1, 1))
print(sobol_seq.i4_sobol_generate(1, 1, 2))
print(sobol_seq.i4_sobol_generate(1, 1, 3))

produced output

[[0.5]]
[[0.5]]
[[0.5]]
[[0.5]]

Code from http://people.sc.fsu.edu/~jburkardt/py_src/sobol/sobol.html, sobol_lib.py behaves reasonable (well, except first point).

Well, enclosed code looks like it might work, keeping seed together with sampled array. Slow, though...

import sobol_seq
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

def i4_sobol_generate_std_normal(dim_num, seed, size=None):
    """
    Generates multivariate standard normal quasi-random variables.
    Parameters:
      Input, integer dim_num, the spatial dimension.
      Input, integer n, the number of points to generate.
      Input, integer seed, initial seed
      Output, real np array of shape (n, dim_num).
    """

    if size is None:
        q, seed = sobol_seq.i4_sobol(dim_num, seed)
        normals = norm.ppf(q)
        return (normals, seed)

    if isinstance(size, int) or isinstance(size, np.int32) or isinstance(size, np.int64) or isinstance(size, np.int16):
        rc = np.empty((dim_num, size))

        for k in range(size):
            q, seed = sobol_seq.i4_sobol(dim_num, seed)
            rc[:,k] = norm.ppf(q)

        return (rc, seed)
    else:
        raise ValueError("Size type is not recognized")

    return None


seed = 1
x, seed = i4_sobol_generate_std_normal(1, seed)
print(x)

x, seed = i4_sobol_generate_std_normal(1, seed)
print(x)

seed = 1
x, seed = i4_sobol_generate_std_normal(1, seed, size=10)
print(x)

x, seed = i4_sobol_generate_std_normal(1, seed, size=1000)
print(x)

hist, bins = np.histogram(x, bins=20, range=(-2.5, 2.5), density=True)
plt.bar(bins[:-1], hist, width = 0.22, align='edge')

plt.show()

Here is the picture

这篇关于Python 中的 Monte Carlo 模拟使用使用 sobol 序列的准随机标准正态数给出了错误的值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！