大学计算-第8章-随机数与采样-课后练习II-随机游走

第1关：带漂移的一维随机游走

import random
random.seed(10)

def random_walk1D_drift(np, ns, r):
    #---------begin-----------
    pos = [0] * np
    for s in range(ns):
        for p in range(np):
            if random.random() < r:
                pos[p] += 1
            else:
                pos[p] -= 1
    return sum(pos)/np

第2关：带漂移的一维随机游走-向量版

import numpy
numpy.random.seed(10)


def random_walk1D_drift(np, ns, r):

    positions = numpy.zeros(np)     # 粒子位置初始化（起点为0）
    #---------begin-----------
    for s in range(ns):
        for p in range(np):
            if numpy.random.random() < r:
                positions[p] += 1
            else:
                positions[p] -= 1
    #---------end-----------

    return numpy.mean(positions)

第3关：有墙的二维随机游走

import numpy 
numpy.random.seed(10)

# 2维有墙随机游走
def random_walk2D_barrier(np, ns):
    # 初始化粒子位置
    X = numpy.zeros(np)
    Y = numpy.zeros(np)
    # 定义矩形区域 A
    xL, yL, xH, yH = -numpy.sqrt(ns), -numpy.sqrt(ns), numpy.sqrt(ns), numpy.sqrt(ns)
    for s in range(ns):
        for p in range(np):
        # 产生方向
            d = numpy.random.randint(1,5)
            s = 1
            # 根据方向和步长移动粒子
            if d==1: #北
                Y[p] += 1
                Y[p] = min(Y[p],yH)
            elif d==2: #南
                Y[p] -= 1
                Y[p] = max(Y[p],yL)
            elif d==3: #西
                X[p] -= 1
                X[p] = max(X[p],xL)
            else: #东
                X[p] += 1
                X[p] = min(X[p],xH)
    # 计算平均位置并返回
    mx, my = X.mean(), Y.mean()
    return mx, my

np = 1000             # 粒子数量

for ns in [100, 200]:   # 步数

    mx, my = random_walk2D_barrier(np, ns)

    #  打印结果

    print('粒子随机游走%d步后平均位置: (%.5f, %.5f)' %(ns, mx, my))

第4关：模拟股票价格

import numpy as np

np.random.seed(10)
import matplotlib.pyplot as plt


def simulate(p0, mu, sigma, T, N):
    dt = T / N
    prices = np.zeros(N + 1)
    prices[0] = p0

    for i in range(1, N + 1):
        drift = mu * prices[i - 1] * dt
        shock = sigma * prices[i - 1] * np.random.normal(scale=np.sqrt(dt))
        prices[i] = prices[i - 1] + drift + shock

    return prices


def draw_picture(prices):
    plt.figure(figsize=(8, 4))
    plt.plot(prices)
    plt.axis([-100, 5200, 8, 29])
    plt.savefig('src/step4/student/result.png')


N = 5000
T = 180
x0 = 10

prices = simulate(x0, 0.01, 0.03, T, N)
draw_picture(prices)