views 1743 words

Intro

When you have a range of values as a result, you are beginning to understand the risk and uncertainty in the model. The key feature of a Monte Carlo simulation is that it can tell you – based on how you create the ranges of estimates – how likely the resulting outcomes are.

思想

1. 所求解的问题本身具有内在的随机性，借助计算机的运算能力可以直接模拟这种随机的过程
1. E.g. 求π
2. 所求解问题可以转化为某种随机分布的特征数
1. E.g.通过随机抽样的方法，以随机事件出现的频率估计其概率，
2. 或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解

Monte Carlo methods vary, but tend to follow a particular pattern:

1. Define a domain of possible inputs.
2. Generate inputs randomly from a probability distribution over the domain.
3. Perform a deterministic computation on the inputs.
4. Aggregate the results.

For example, consider a circle inscribed in a unit square. Given that the circle and the square have a ratio of areas that isπ/4, the value of π can be approximated using a Monte Carlo method:

1. Draw a square, then inscribe a circle within it.
2. Uniformly scatter objects of uniform size over the square.
3. Count the number of objects inside the circle and the total number of objects.
4. The ratio of the two counts is an estimate of the ratio of the two areas, which is π/4. Multiply the result by 4 to estimate π.

import random
from math import sqrt, pi

x = y = inn = out = 0.0

tolerance = 1e-5
piApproximation = 0.0
# for i in range(30000):
while abs(pi - piApproximation) > tolerance:
x = random.random()
y = random.random()
# print(x, y)
# print(sqrt(x * x + y * y))
# 假设正方形边长为2.0
# 点到圆心的距离需要小于1.0, 才证明点在圆内
if sqrt(x**2 + y**2)<1.0:
# if (x**2 + y**2) <= 1.0:
inn += 1.0
else:
out += 1.0
piApproximation = 4.0 * (inn / (inn+out))

print(inn)
print(out)
print(piApproximation)
# print(inn/30000 * 4)
print(inn / (out + inn) * 4)

π的求解

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle

# 随机值数量 半径 圆心位置
# n = 100000
n = 10
r = 1.0
a,b, = 0.0, 0.0

np.random.seed(42)

y_max = b + r  # 1
y_min = b - r  # -1
x_max = a + r  # 1
x_min = a - r  # -1

# np.random.uniform()函数生成均匀分布的随机数
x = np.random.uniform(x_min,x_max,n) # [-0.25091976  0.90142861  0.46398788  0.19731697, ...
y = np.random.uniform(y_min,y_max,n) # [-0.95883101  0.9398197   0.66488528 -0.57532178, ....
print(x)
print('======')
print(y)
print('=======')
distance = np.sqrt((x-a)**2+(y-b)**2) # 点到圆心的距离
print(distance) # [0.99111939 1.30224215 0.81077567 0.608218, ...

count = sum(np.where(distance <= r,1,0))
print(count)

pai = 4 * count / n

print('π的值为:',pai)

fig = plt.figure(figsize = (12,12))
plt.scatter(x,y,marker='.',s = 10,alpha = 0.5)
circle = Circle(xy = (a,b),radius = r, alpha = 0.5 ,color = 'gray')
plt.grid()
yuan.grid(color='k' ,linestyle='--', alpha = 0.2,linewidth=2) # 设置背景网格线的样式
# 绘制圆形

C# code

using System;

class MainClass {
public static void Main (string[] args) {
// Construct a random number generator that generates random number greater than or equal to 0.0, and less than 1.0
Random rand = new Random();

// Approximate pi to within 5 digits.
double tolerance = 1e-5;

double piApproximation = 0;
int total = 0;
int numInCircle = 0;
double x, y; // Coordinates of the random point.

// Generate random points until our approximation within the desired tolerance.
while ( Math.Abs( Math.PI - piApproximation ) > tolerance )
{
x = rand.NextDouble();
y = rand.NextDouble();
if ( x * x + y * y <= 1.0 ) // Is the point in the circle?
{
++numInCircle;
}
++total;
piApproximation = 4.0 * ( (double) numInCircle / (double) total );
}

Console.WriteLine();

Console.WriteLine( "Pi calculated to within {0} digits with {1} random points.",
-Math.Log10( tolerance ), total );
Console.WriteLine( "Approximated Pi = {0}", piApproximation );
}
}


求定积分

通用情况

1. 产生服从分布律FXFX的随机变量Xi(i=1,2,⋯,N)Xi(i=1,2,⋯,N)
2. 计算均值

$$\overline I = \frac{1}{N} \sum_{i = 1}^N g^*(X_i)$$ 此时有 $$\overline I \approx I$$

代码

n = 100000
x_min, x_max = 0.0, 1.0
y_min, y_max = 0.0, 1.0

count = 0
for i in range(0, n):
x = random.uniform(x_min, x_max)
y = random.uniform(y_min, y_max)
# x*x > y，表示该点位于曲线的下面。所求的积分值即为曲线下方的面积与正方形面积的比。
if x*x > y:
count += 1

integral_value = count / float(n)
print('定积分的结果为: ', integral_value)

fig = plt.figure(figsize = (12,12))
axes.plot(x, y,'o',markersize = 1)
# 绘制散点图

xi = np.linspace(0,1,100)
yi = xi ** 2
plt.plot(xi,yi,'--r')
plt.fill_between(xi, yi, 0, color ='gray',alpha=0.5,label='area')
plt.grid()
# 绘制 y = x**2 面积图

证券市场

• 如果交易冷清，你会以平均价11元，卖出5万股。
• 如果交易活跃，你会以平均价8元，卖出10万股。
• 如果交易温和，你会以平均价10元，卖出7.5万股。

As a result, you expect to sell 75,000 units (i.e., (50,000+75,000+100,000)/3 = 75,000) at an average selling price of \$9.67 per unit (i.e., (\$11+\$10+\$8)/3 = $9.67) calculates Net Profit based on average sales volume, average selling price, and average unit cost. Flawed Average = 75,000 * (\$9.67 - \$6.5) - 12,000 = 117,500 (The true average Net Profit is closer to$93,000, the average will nearly always be much less than \$117,750) 已知你的成本(unit cost)在每股5.5元到7.5元之间，平均是6.5元。请问接下来的交易，你的净利润会是多少？ 取1000个随机样本，每个样本有两个数值：一个是证券的成本（5.5元到7.5元之间的均匀分布），另一个是当前市场状态（冷清、活跃、温和，各有三分之一可能 • 1000次抽样，每次按均匀分布随机选定一个成本值、同时按1/3的概率选中一个市场情境，然后拿这两个参数可计算出净利润、平均利润，千次加总并取均值就是结果。1000次抽样的技术本身就是蒙特卡罗模拟 Net profit will be calculated as Net Profit = Sales Volume* (Selling Price - Unit cost) - Fixed costs. Fixed costs (for overhead, advertising, etc.) are known to be$120,000

### way1
avg_cost = 6.5
sum_profit = 0
for i in range(1000):
# cost = random.uniform(5.5,7.5)
avg_sell = random.choice([(11,50000),(8,100000),(10,75000)])
avg_profit = (avg_sell[0] - avg_cost) * avg_sell[1] - 120000  # 120000 is fixed costs
sum_profit += avg_profit

avg_profit = sum_profit / 1000
print(avg_profit)  # 92775.0

### way2
net_profit = 0
for i in range(1000):
cost = random.uniform(5.5,7.5)
profit = (11 - cost)*50000 + (8 - cost)*100000 + (10 - cost)*75000 - 120000*3  # 120000 is fixed costs
net_profit += profit
res = net_profit / (3*1000)
print(res)   # 92025.0

C#:

using System;

class MainClass {
public static void Main (string[] args) {

Random rand = new Random();
double net_profit = 0.0;

for(int i =0; i< 1000; i++){
double cost = (rand.NextDouble() * 2.0) + 5.5;
// Console.WriteLine(cost);
double profit = (11 - cost)*50000 + (8 - cost)*100000 + (10 - cost)*75000 - 120000*3;

net_profit += profit;

}
net_profit = net_profit / (3*1000);
Console.WriteLine( "Net profit = {0}", net_profit );

}
}

using System;

class MainClass {
public static void Main (string[] args) {

Random rand = new Random();
double net_profit = 0.0;
double avg_cost = 6.5;
int[][] pos = new int[3][]{new int[]{11,50000}, new int[]{8,100000}, new int[]{10,75000}};

// int temp = rand.Next(pos.Length);
// Console.WriteLine(pos);
// Console.WriteLine(pos.Length);
// Console.WriteLine(temp);
// Console.WriteLine(pos[temp][0]);

for(int i =0; i< 1000; i++){
int temp = rand.Next(pos.Length);
double cost = pos[temp][0];
double vol = pos[temp][1];
double profit = (cost - avg_cost) * vol - 120000;

net_profit += profit;

}
net_profit = net_profit / 1000;
Console.WriteLine( "Net profit = {0}", net_profit );

}
}


stock price prediction using monte carlo

1. stock price data in dataframe
2. calculate log return
• log_return = log(1 + data.pct_change())
• pct_change() stands for percentage change <==> moving average
3. calculate mean, variance and std
• u = log_return.mean()
• var = log_return.var()
• stdev = log_return.std()
4. compute the drift component (the best approximation of future rates of return of the stock)
• $drift = u - \frac{1}{2} \cdot var$
• average log return - half its variance
5. Brownian motion
• first component:
• $r = drift + stdev * e^r$
• Brownian motion comprises the sum of the drift and standard deviation of adjusted by e to the power of r
• second component:
• a random variable Z: Z corresponds to the distance between the mean and the events, expressed as the number of standard deviations
• E.g. norm.ppf(0.95) = 1.65
• if an event has a 95% chance of occurring, the distance between this event and mean will approximately 1.65 standard deviations
• Z = spicy.stat.norm.ppf(np.random.rand(10,2))
• use the probabilities generated by the rand and convert them into distances from the mean zero as measured by the number of standard deviations
• $daily_returns = e^r$
• $r = drift + stdev * z$
• e is Euler’s number / natural number
• Next Day’s Price = Today’s Price x e^(Drift+Random Value)
• time_interval is 1000 (since it is forecasting the stock price for the upcoming 1000days)
• iteration is 10 (ask computer to produce 10 series of future stock price prediction)
• daily_returns = np.exp(drift + stdev * norm.ppf(np.random.rand(time_interval, iteration)))
• 1000rows x 10cols
• S0 = data.iloc[-1] which is the last adjusted closing price of the stock
• price_list= np.zeros_like(daily_returns) – create a price_list, the same dimension as the daily_returns matrix with all 0.0
• price_list[0] = S0 - Set the values on the first row of the price_list array equal to S0
• Create a loop in the range (1 to t_intervals) that reassigns to the price in time t the product of the price in day (t-1) with the value of the daily returns in t
• for t in range(1, t_intervals): price_list[t] = price_list[t – 1] * daily_returns[t]

code:

#Import Libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
%matplotlib inline
﻿
#Settings for Monte Carlo asset data, how long, and how many forecasts
ticker = 'XTZ_USD' # ticker
t_intervals = 30 # time steps forecasted into future
iterations = 25 # amount of simulations
#Acquiring data
data = pd.read_csv('XTZ_USD Huobi Historical Data.csv',index_col=0,usecols=['Date', 'Price'])
data = data.rename(columns={"Price": ticker})
#Preparing log returns from data
log_returns = np.log(1 + data.pct_change())
#Plot of asset historical closing price
data.plot(figsize=(10, 6));

#Plot of log returns
log_returns.plot(figsize = (10, 6))

#Setting up drift and random component in relation to asset data
u = log_returns.mean()
var = log_returns.var()
drift = u - (0.5 * var)
stdev = log_returns.std()
daily_returns = np.exp(drift.values + stdev.values * norm.ppf(np.random.rand(t_intervals, iterations)))
#Takes last data point as startpoint point for simulation
S0 = data.iloc[-1]
price_list = np.zeros_like(daily_returns)
price_list[0] = S0
#Applies Monte Carlo simulation in asset
for t in range(1, t_intervals):
price_list[t] = price_list[t - 1] * daily_returns[t]
﻿
﻿﻿#Plot simulations
plt.figure(figsize=(10,6))
plt.plot(price_list);

Ref

https://www.shuzhiduo.com/A/VGzlbNGV5b/

https://www.shuzhiduo.com/A/x9J2QZ1Vz6/

马尔可夫链蒙特卡洛采样 (Markov Chain Monte Carlo - MCMC)

MCMC 方法用于通过在概率空间中进行随机采样以近似地得出某一感兴趣参数的后验分布

https://baijiahao.baidu.com/s?id=1628127197886312445&wfr=spider&for=pc

https://baijiahao.baidu.com/s?id=1588750123445231830&wfr=spider&for=pc