import numpy as np
import matplotlib.pyplot as plt

# 设置参数
T = 1.0      # 终止时间
N = 1000     # 时间步数
dt = T / N   # 每一步的时间间隔
sigma = 0.2  # 波动率
S0 = 100     # 初始价格
A = 1.0      # 基准强度
k = 0.5      # 衰减系数
delta = 1.0  # 价格距离

# 生成布朗运动路径
t = np.linspace(0, T, N+1)  # 时间点
W = np.random.randn(N)      # 生成N个标准正态分布的随机数
W = np.insert(W, 0, 0)      # 将初始值0插入到W的开头
W = np.cumsum(W) * np.sqrt(dt)  # 计算布朗运动的累加和并乘以sqrt(dt)

# 计算参考价格路径
S = S0 + sigma * W

# 计算捕获流强度
Lambda_plus = A * np.exp(-k * delta + k * (S - S[0]))
Lambda_minus = A * np.exp(-k * delta - k * (S - S[0]))

# 绘制参考价格和捕获流强度
fig, ax1 = plt.subplots(figsize=(10, 6))

# 绘制参考价格
ax1.set_xlabel('Time')
ax1.set_ylabel('Reference Price', color='tab:blue')
ax1.plot(t, S, label='Reference Price', color='tab:blue')
ax1.tick_params(axis='y', labelcolor='tab:blue')

# 创建第二个y轴
ax2 = ax1.twinx()
ax2.set_ylabel('Intensity', color='tab:red')
ax2.plot(t, Lambda_plus, label=r'$\Lambda^+(\delta, t, t + \Delta T)$', color='tab:red')
ax2.plot(t, Lambda_minus, label=r'$\Lambda^-(\delta, t, t + \Delta T)$', color='tab:green')
ax2.tick_params(axis='y', labelcolor='tab:red')

# 设置图例
fig.tight_layout()
plt.title('Reference Price and Captured Flow Intensities')
fig.legend(loc='upper left', bbox_to_anchor=(0.1,0.9))
plt.grid(True)
plt.show()