R 到 Python 随机过程转换

Posted 2023-02-14

【中文标题】R 到 Python 随机过程转换

【英文标题】：R to Python stochastic process translation

【发布时间】：2022-01-17 18:48:22

【问题描述】：

我想把下面的 R 代码翻译成 Python。

这主要是一个随机过程，我需要将其翻译成 Python。

代码实现了具有两个波动阶段的跳跃过程的马尔可夫链模拟。

set.seed(42)

nSim   <- 1E5
tau    <- 3
K      <- 105
S0     <- 100
rf     <- 0.05
vol_lo <- 0.25
vol_hi <- 0.75
lambda <- c(3,2) # away-from-lo, away-from-hi

sim_time_in_lo <- function(state0)
  t <- 0
  s <- state0
  time_lo <- 0
  while (t<tau)
    dt <- rexp(n=1,lambda[s])
    if ((t+dt)>tau) dt <- tau - t
    if (s==1)time_lo <- time_lo + dt 
    if (s==1)s<-2 else s <-1
    t <- t+ dt
  
  time_lo

tau_lo     <- sapply(1:nSim,function(i)sim_time_in_lo(1))
tau_hi     <- tau - tau_lo
total_var  <- tau_lo * vol_lo^2 + tau_hi * vol_hi^2
drift      <- rf * tau-0.5*total_var
randomness <- sqrt(total_var)*rnorm(nSim,)

Python 尝试如下：

我认为我在 pandas 中使用了类似 sapply 的功能。

import numpy as np
import pandas as pd


np.random.seed(42)

nSim   = 1000
tau    = 3
K      = 105
S0     = 100
r      = 0.05
vol_lo = 0.25
vol_hi = 0.75
lambd = [3,2] # away-from-lo, away-from-hi

def sim_time_in_lo(state0):
    t = 0
    s = state0
    time_lo = 0
    while (t<tau):
        dt = np.random.exponential(scale = lambd[s],size = nSim)
        if (t+dt)>tau:
            dt = tau - t
        elif s == 1:
                time_lo = time_lo + dt
        elif s == 1:
            s = 2
        else : s = 1
        t = t+ dt
        time_lo

df = pd.DataFrame([sim_time_in_lo(1) for i in range(1,nSim)]);df
tau_hi     = tau - tau_lo
total_var  = tau_lo * vol_lo**2 + tau_hi * vol_hi**2
drift      = r * tau-0.5*total_var
randomness = np.sqrt(total_var)*np.random.normal(size=nSim)

但我收到一个错误

df = pd.DataFrame([sim_time_in_lo(1) for i in range(1,nSim)]);df

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

我的错误是什么？

【问题讨论】：

您的列表理解如何？你能给我们提供：[sim_time_in_lo(1) for i in range(1,nSim)][0:5]

我不太擅长 R，但我可以告诉你 Python 的错误是什么——你正在尝试将向量 (t+dt) 与标量 (tau) 进行比较。问哪个更大不是一个格式正确的问题。我需要更好地理解 R 代码在做什么

@MathiasByskov 它也没有

【参考方案1】：

以下是一些使其运行的改进：

nSim = int(1e5)
tau = 3
K = 105
S0 = 100
r = 0.05
vol_lo = 0.25
vol_hi = 0.75
lambd = [3, 2]  # away-from-lo, away-from-hi


def sim_time_in_lo(state0):
    t = 0
    s = state0
    time_lo = 0
    while t < tau:
        dt = np.random.exponential(size=1, scale=1/lambd[s])
        if (t + dt) > tau:
            dt = tau - t
        if s == 0:
            time_lo = time_lo + dt
            s = 1
        else:
            s = 0
        t = t + dt
    if isinstance(time_lo, np.ndarray):
        return time_lo[0]
    return time_lo


tau_lo = np.array([sim_time_in_lo(0) for i in range(nSim)])
tau_hi = tau - tau_lo
total_var = tau_lo * vol_lo ** 2 + tau_hi * vol_hi ** 2
drift = r * tau - 0.5 * total_var
randomness = np.sqrt(total_var) * np.random.normal(size=nSim)

nSim 是 1e5，而不是 1000。这种语法也适用于 Python，但您需要转换为 int。 与 R 不同，Python 使用零索引，因此状态需要为 0 和 1 才能从 lambd 获取值。 指数的 size 参数应为 1，以返回标量，而不是 nSim - 这就是您收到错误的原因。 您必须使用显式的 return 语句，这与 R 不同（否则 Python 函数会返回 None）。 修复了一些奇怪的缩进，我认为这些缩进也会导致问题。 首先将time_lo 初始化为大小为1 的数组，否则有时会返回标量。然后在最后从 size-1 数组中提取标量。 使用数组而不是数据框。 （现在不需要导入 pandas。） 注意 np.random.exponential的scale参数是$\beta$，即$\frac1\lambda$。 R 的rexp 的rate 是$\lambda$。

【讨论】：

我收到了错误TypeError: 'int' object is not subscriptable

错误是指tau_lo = np.array([sim_time_in_lo(0)[0] for i in range(nSim)]) 行

好的，我编辑了return语句，再试一次

【参考方案2】：

您的错误在于if (t+dt)>tau。 dt 是一个数组，它使 (t+dt)>tau 成为一个布尔值数组。

使用((t+dt)>tau).any() 或((t+dt)>tau).all() 赋予if 语句含义。 .all() 表示您希望 t+dt 的每个值都大于 tau，.any() 表示如果只有一个值更大，那就足够了。

【讨论】：

我将其更改为.any()，但收到错误TypeError: loop of ufunc does not support argument 0 of type float which has no callable sqrt method