MATLAB 笔记：蒙特卡洛方法

由polll创建，最终由polll更新于2024-05-20 03:18 被浏览 30 用户

蒙特卡洛(Monte Carlo)方法是一种基于随机数的计算方法。这一方法源于美国在二战期间研制原子弹的“曼哈顿计划”，该计划的主持人冯诺依曼用摩纳哥驰名世界的赌城Monte Carlo来命名这个方法，因此称之为Monte Carlo方法。

Monty Hall Problem，也称为三门问题，是一个源自博弈论的数字游戏问题。这个游戏的玩法是：参赛者面前有三扇关闭的门，其中一扇门的后面藏有一辆汽车，而另外两扇门的后面各自藏有一只山羊。参赛者从三扇门中随机选取一扇，若选中后面有车的门就可以赢得汽车。当参赛者选定了一扇门，但尚未开启它的时候，节目主持人会从剩下两扇门中打开一扇藏有山羊的门，然后问参赛者要不要更换自己的选择，选取另一扇仍然关上的门。这个游戏涉及到的问题是，参赛者更换自己的选择是否会增加赢得汽车的概率？


Vos Savant给出的解答是参赛者应该更换选择。根据全概率公式，如果参赛者坚持自己原来的选择，那么赢得车的概率是1/3。而如果更换选择，概率变为2/3。


我们还可以通过蒙特卡洛方法求解参赛者更换选择之后赢得汽车的概率。设两只羊的编号分别为1和2，汽车编号为3.现在从数字1，2，3中随机选取一个数字，若一开始选中1或2，则更换选择后选中3，即赢得汽车；若一开始选中3，则更换选择后选中1获，即得不到汽车。将这样的实验重复进行n次，记录一开始选中1或2的次数m，从而可以确定更换选择后赢得汽车的频率m/n。由大数定律可知当实验次数n增大时，频率m/n趋近于更换选择后赢得汽车的概率。

function p=SheepAndCar(n)
% p=SheepAndCar(n)，利用蒙特卡洛方法求解monty hall问题，
%求参赛者更换选择之后赢得汽车的概率p。
for i=1:length(n)
    x=randsample(3,n(i),'true'); %随机抽样
    p(i)=sum(x~=3)/n(i);
end

p=SheepAndCar([10,100,1000,10000,100000,1000000]) %求概率模拟值向量
p =

    0.2000    0.6900    0.6710    0.6722    0.6668    0.6664

由以上结果可以看到，随着随机抽样次数的增大，所求概率的模拟值逐渐趋近于理论值2/3。

用蒙特卡洛方法求积分


用随机投点法求解。如图，由两条虚线(x=1,y=1)和坐标轴围城了一个边长为1的正方形。在这个正方形内随机投点。所投点落到阴影区域内的概率等于A面积与正方形面积之比。当随机投点总数足够大时，用落到阴影区域内的点的频率m/n近似值。matlab程序如下：

function[S0,Sm]=quad1mont(n)
%[S0,Sm]=quad1mont(n),求曲线y=sqrt(x)与直线y=x所围城的阴影区域的面积理论值S0
%与monte carlo理论值sm。
%S0=int('sqrt(x)-x',0,1); %面积的理论值（解析解）
S0=quad(@(x)sqrt(x)-x,0,1)%面积的理论值（数值解）
%计算阴影区域面积的蒙特卡洛模拟值
for i=1:length(n)
x=rand(n(i),1); %点的横坐标
y=rand(n(i),1); %点的纵坐标
m=sum(sqrt(x)>=y&y>=x); %落到阴影区域内点的频数
Sm(i)=m/n(i);%落到阴影区域内点的频率，即概率的模拟值
end
%执行程序
[v0,vm]=quad1mont([10,100,1000,10000,100000,1000000])


二重积分

求球体 被圆柱面 所截得的（含在圆柱面内的部分）立体的体积。


function[v0,vm]=quad2mont(n)
%v0=32*(pi/2-2/3)/3; %体积的理论值
v0=4*quad2d(@(x,y)sqrt(4-x.^2-y.^2),0,2,0,@(x)sqrt(1-(1-x).^2));
%求体积的蒙特卡洛模拟值
for=i=1:length(n)
 x=2*rand(n(i),1);    %点的x坐标
y=rand(n(i),1); %点的y坐标
z=2*rand(n(i),1);%点的z坐标
%落到区域T内的点的频数
m=sum((x.^2+y.^2+z.^2<=4)&((x-1).^2+y.^2<=1));
vm(i)=16*m/n(i); %落到所求立体内的点的频率
end

模拟股票价格。初始价格为20，股票价格变化根据集合布朗运动过程计算：


mu为期望收益率，sigma为波动项， 是一个几何布朗运动。进行10000次模拟，每次计算90天内价格的变化，之后将每天每次实验得出的价格求平均值，求得经过10000次模拟经验总结的股票价格变化以及价格分布。在信用评级、风险管理领域多有应用。

clear
dt=1/365.0;                   % 一天的年单位时间
S0=20;                       % 股票在初始时刻的价格，程序中假设
r=0.031;                      % 期望收益率
sigma=0.6;                    % 波动率=0.6
expTerm=r*dt;                 % 漂移项dt
stddev=sigma*sqrt(dt);       % 波动项o:dz(t) 
nDays1=90;                    % 要模拟的总天数
for nDays=1:nDays1            % nDays表示时刻t
nTrials=10000;                % 模拟次数
for j=1:nTrials 
n = randn(1,nDays);           %生成nDays个标准正态分布随机数
S=S0; 
for i=1:nDays 
dS = S*(expTerm+stddev*n(i));   % 模拟计算股票价格的增量
S=S+dS;                      %计算股票价格
end
S1(nDays,j)=S;               % 将每天的股票模拟价格数据记录在S1中
end
end
S2=mean(S1');                 % 计算每天模拟的股票价格的均值，作为价格的估值 
plot(S2','-o')                    % 90天期间股票价格估值的曲线图
figure(2)
hist(S1(90,:),0:0.5:65)      %第90天的股票价格模拟的直方图



用Black-Scholes期权定价公式模拟股票价格


T=1;
S0=48; % underlying asset（股票为主）期初价格
r=0.05; % risk free rate 无风险利率
NStep=200; % number of steps 将时间分成200份
deltaT=T/NStep; % Time Step 
Sigma=0.1; % Strand Derivation 
Rep=20000; % Number of replication 重复试验次数
mu=( r-Sigma^2/2 )*deltaT;
sigmadelta=Sigma*sqrt(deltaT);
% 计算增量
Increments=mu+sigmadelta.*randn(NStep,Rep); 增量为(r-sigma^2/2)*(T-t)加上标准差*
%0，1正态分布的数，一共200个20000次模拟
% 计算对数价格
LogPrice=cumsum([log(S0)*ones(1,Rep); Increments]); 累加增量之后返回1到20000次的价格
PricePath=exp(LogPrice);

%对价格路径画图
plot(PricePath)
title('Simulation of Price Path')
ylabel('Stock Price');
xlabel('Time Step');
xlim([0,NStep]);

重复20次试验之后的结果


重复20000次试验之后的结果