如何在Python中为布朗运动蒙特卡洛模拟加入自相关修正？

DevOps 小满

2026-7-12

问题描述

本人是Python及论坛新手，若有理解偏差请见谅。现发现连续值间存在负一阶自相关，希望将其纳入布朗运动蒙特卡洛（BMMC）模拟，认为此举可缩小随时间变化的置信区间。以下为无自相关的模拟代码，参考了相关论文但未能复现结果，现提出两点疑问：

上述理解是否正确，该方法是否适用？若不适用，有哪些可纳入负一阶自相关的替代方案？
如何在现有代码中正确实现最优方法以复现论文类似结果？

# Data
data = pd.DataFrame([3.1, 4.2, 3.7, 5.1, 4.3, 6.6, 8.4, 8.1, 7.5, 11.6, 10.2, 13.2])

# Calculate the logarithmic, differenciated y of our data.
df = np.log(data).diff().dropna()

# Extract the last available value 
S0 = data.iloc[len(df)]


mu = df.mean()
var = df.var()
desvest = df.std()
n_simulations = 2000
time_units = np.arange(9)
conf_level = 0.8
z_score = NormalDist().inv_cdf((1 + conf_level) / 2.)


# Determine the drift (slope)
drift = mu - (0.5*var)

# Determine the epsilon
epsilon = norm.ppf(np.random.rand(len(time_units), n_simulations))

# Determine y
y =  drift.values + desvest.values * epsilon

# Prepare the series
y_interval_1 = np.zeros(len(time_units))
y_interval_2 = np.zeros(len(time_units))
expected_y = np.zeros(len(time_units))

# Compute the confidence interval logic
for t in range(1, len(time_units)):
    y_interval_1[t] = drift*(t-time_units[0]) + desvest *z_score*np.sqrt(t - time_units[0])
    y_interval_2[t] = drift*(t-time_units[0]) + desvest *(-1*z_score)*np.sqrt(t - time_units[0])

# Prepare the outcome series   
S = np.zeros_like(y)
S_interval_1 = np.zeros_like(y_interval_1)
S_interval_2 = np.zeros_like(y_interval_2)

# Set initial value(s) - last available value of inserted series, named S0
S_interval_1[0] = S0
S_interval_2[0] = S0
S[0] = S0
expected_y[0] = S0

# Simulate the individials trails by using a for loop
for t in range(1, len(time_units)):
    S[t] = S[t-1]*np.exp(y[t])
    S_interval_1[t] = S0*np.exp(y_interval_1[t])
    S_interval_2[t] = S0*np.exp(y_interval_2[t])
    expected_y[t] = expected_y[t-1]*np.exp(mu.values)

解答

1. 理解合理性与替代方案

理解正确性

你的核心逻辑是对的：负一阶自相关意味着相邻收益会反向抵消，长期来看收益的方差累积速度会慢于无自相关的情况，因此对应的置信区间会更窄。但标准几何布朗运动（GBM）假设收益独立同分布，无法直接纳入自相关结构，所以直接修改现有GBM代码无法实现你的需求。

替代方案

针对一阶负自相关的场景，推荐以下几种模型：

AR(1)模型（一阶自回归）：直接对对数收益建模，捕捉一阶自相关关系，是最贴合你需求的轻量方案。
ARMA-GARCH模型：如果除了均值的自相关，波动率也存在聚类特征，可同时建模均值自相关和波动率异方差。
带自相关扰动的随机游走：修改GBM的扰动项，让其服从AR(1)过程，保留GBM的框架但加入自相关约束。

2. 代码实现（基于AR(1)模型的最优方案）

我们选择AR(1)模型来修改现有代码，步骤分为：估计AR(1)参数、生成带自相关的模拟收益、重新计算置信区间。

修改后的完整代码

import pandas as pd
import numpy as np
from scipy.stats import norm, NormalDist
from statsmodels.tsa.arima.model import ARIMA

# 1. 数据预处理
data = pd.DataFrame([3.1, 4.2, 3.7, 5.1, 4.3, 6.6, 8.4, 8.1, 7.5, 11.6, 10.2, 13.2])
log_returns = np.log(data).diff().dropna().values.flatten()  # 转为一维数组方便建模
S0 = data.iloc[-1].values[0]  # 取最后一个观测值

# 2. 估计AR(1)模型参数
ar1_model = ARIMA(log_returns, order=(1, 0, 0)).fit()
mu_ar1 = ar1_model.params[0]  # 均值项
rho_ar1 = ar1_model.params[1]  # 一阶自相关系数
resid_var = ar1_model.resid.var()  # 残差方差
resid_std = np.sqrt(resid_var)

# 3. 模拟参数设置
n_simulations = 2000
time_units = np.arange(9)
conf_level = 0.8
z_score = NormalDist().inv_cdf((1 + conf_level) / 2.)

# 4. 生成带自相关的模拟收益序列
# 初始化收益矩阵，每行对应一个时间步，每列对应一个模拟路径
sim_returns = np.zeros((len(time_units), n_simulations))
# 第一个时间步的收益：从AR(1)的平稳分布生成
sim_returns[0] = np.random.normal(loc=mu_ar1 / (1 - rho_ar1), scale=resid_std / np.sqrt(1 - rho_ar1**2), size=n_simulations)
# 后续时间步按照AR(1)生成：r_t = mu + rho*(r_{t-1} - mu) + eta_t
for t in range(1, len(time_units)):
    sim_returns[t] = mu_ar1 + rho_ar1 * (sim_returns[t-1] - mu_ar1) + np.random.normal(loc=0, scale=resid_std, size=n_simulations)

# 5. 模拟价格路径
S = np.zeros_like(sim_returns)
S[0] = S0
for t in range(1, len(time_units)):
    S[t] = S[t-1] * np.exp(sim_returns[t])

# 6. 计算带自相关的置信区间
# 长期方差计算公式：Var(r_1+...+r_T) = T*sigma_eta²/(1-rho²) + 2*rho*sigma_eta²/(1-rho²)*(T - (1-rho^T)/(1-rho))
expected_log_price = np.log(S0) + mu_ar1 / (1 - rho_ar1) * time_units
var_total = np.zeros(len(time_units))
for t in range(len(time_units)):
    if t == 0:
        var_total[t] = 0
    else:
        T = t
        term1 = T * resid_var / (1 - rho_ar1**2)
        term2 = 2 * rho_ar1 * resid_var / (1 - rho_ar1**2) * (T - (1 - rho_ar1**T)/(1 - rho_ar1))
        var_total[t] = term1 + term2
std_total = np.sqrt(var_total)

S_interval_1 = np.exp(expected_log_price + z_score * std_total)
S_interval_2 = np.exp(expected_log_price - z_score * std_total)
expected_y = np.exp(expected_log_price)

关键修改说明

AR(1)参数估计：用statsmodels的ARIMA模型（order=(1,0,0)即AR(1)）估计均值、自相关系数和残差方差。
收益生成逻辑：不再生成独立正态扰动，而是按照AR(1)的递推公式生成带负自相关的收益序列。
置信区间计算：考虑自相关后的方差累积公式，替代原GBM的简单t*sigma²，确保置信区间能正确反映负自相关带来的方差缩减。

内容的提问来源于stack exchange，提问作者user22465897

如何在Python中为布朗运动蒙特卡洛模拟加入自相关修正？

1. 理解合理性与替代方案 #

理解正确性 #

替代方案 #

2. 代码实现（基于AR(1)模型的最优方案） #

修改后的完整代码 #