蒙特卡洛方法是一种利用计算机的随机数理论模拟实际的情况的一种方法。今天主要是以实例讲解蒙特卡洛方法的MATLAB编程实现。
实例1
程序
clc;%清除命令行窗口命令
clear all;%清除工作区变量
close all;%关闭图形窗口
L = 1;%正方形长度
n = 100000;%计算随机生成的10万点数 科学计数法
%随机生成[0,1]范围内的随机数
x = unifrnd(0,L,[1 n]);%%unifrnd 连续均匀分布的随机数生成器
y = unifrnd(0,L,[1 n]);
pinshu = sum((x.^2+y.^2<=L));%判断点是否在单位圆内
ix = find(x.^2+y.^2<=L);
%绘图
figure;
plot([0 L L 0],[0 0 L L],'k-');
hold on;
plot(x(ix),y(ix),'r.');
Pai = 4*pinshu/n;
%输出结果
fprintf('蒙特卡洛算法求解圆周率:%.20f\n',Pai);
xlim([0 L*1.2]);
ylim([0 L*1.2]);
legend('正方形区域','单位圆内的点'); 运行结果
实例2
程序
clc;
clear all;
close all;
syms x y
f1 = x.^4+2.*y.^4-81;
f2 = x.^2-y.^2-4;
figure;
ezplot(f1);
hold on;
ezplot(f2);
n = 1e6;
x = unifrnd(-6,6,[1 n]);%unifrnd 连续均匀分布的随机数生成器
y = unifrnd(-6,6,[1 n]);
pinshu = sum((x.^4+2.*y.^4<=81)&(x.^2-y.^2<=4));
area = 12*12*pinshu/n
count = 0;
for i = 1:n
if (x(i)^4+2*y(i)^4)<=81&(x(i)^2-y(i)^2<=4)
count = count +1;
x1(count) = x(i);
y1(count) = y(i);
end
end
hold on;
plot(x1,y1,'r*');
title('面积区域')
%结果area =22.7863 运行结果
实例3
模拟抛硬币，计算正面出现的频率，并画出随着试验次数n的增大，频率和概率的关系图。(unidrnd(N):生成最大值为N的随机正整数。)
程序
clc;
clear all;
close all;
%模拟抛硬币
n = 100:1:1000;
for i = 1:length(n)
num = n(i);
x = rand(1, num);%rand(size) 随机生成0到1之间均匀分布的浮点数
% 假设大于0.5为硬币向上
index = find(x>0.5);
count(i) = length(index);
p(i) = count(i)/num;
end
figure;
plot(n,p,'r-');
xlabel('次数n');
ylabel('正面向上的频率')
grid on;
ylim([0 1]);
disp('随着n不断增大，频率稳定在概率0.5。') 运行结果
实例4
程序
clc;
clear all;
close all;
a = 0;
b = 1;
f =@(x) 4./(1+x.^2);
n = [100 1000 10000 100000]
for i = 1:length(n)
t = rand(1,n(i));%rand(size) 随机生成0到1之间均匀分布的浮点数
x = a+(b-a)*t;
s = sum(f(x));
S(i) = s*(b-a)/n(i);
end
S
error = pi-S 运行结果
n =
100 1000 10000 100000
S =
3.1996 3.1227 3.1372 3.1419
error =
-0.0580 0.0189 0.0043 -0.0003
>> 网络上其他蒙特卡洛计算程序
程序
clear all;
clc;
close all;
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天的股票价格模拟的直方图
运行结果
程序
clc;
clear all;
close all;
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]); 运行结果
本文内容来源于网络，仅供参考学习，如内容、图片有任何版权问题，请联系处理，24小时内删除。
作 者 | 郭志龙
编 辑 | 郭志龙
校 对 | 郭志龙