蒙特卡洛方法是一種利用計算機的随機數理論模拟實際的情況的一種方法。今天主要是以執行個體講解蒙特卡洛方法的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小時内删除。
作 者 | 郭志龍
編 輯 | 郭志龍
校 對 | 郭志龍