《2022年随机过程matlab程序 .pdf》由会员分享，可在线阅读，更多相关《2022年随机过程matlab程序 .pdf（17页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、% PPT 例2 一维正态密度与二维正态密度syms x y; s=1; t=2; mu1=0; mu2=0; sigma1=sqrt(1+s2); sigma2=sqrt(1+t2); x=-6:0.1:6; f1=1/sqrt(2*pi*sigma1)*exp(-(x-mu1).2/(2*sigma12); f2=1/sqrt(2*pi*sigma2)*exp(-(x-mu2).2/(2*sigma22); plot(x,f1,r-,x,f2,k-.) rho=(1+s*t)/(sigma1*sigma2); f=1/(2*pi*sigma1*sigma2*sqrt(1-rho2)*exp
2、(-1/(2*(1-rho2)*(x-mu1)2/sigma12-2*rho*(x-mu1)*(y-mu2)/(sigma1*sigma2)+(y-mu2)2/sigma22); ezsurf(f) -6-4-2024600.050.10.150.20.250.30.35-505-50500.050.10.150.2x44798133900177/281474976710656 exp(-5/2 x2+3 x y-y2)y% % The daily log returns on the stock have a mean of 0.05/year and a standard deviatio
3、n of 0.23/year. These can be converted to rates per trading day by deviding by 253 and sqrt(253), respectively.名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 17 页 - - - - - - - - - Question 1: What is the probability that the value of the stock will be below $9
4、50,000 at the close day of at least one of the next 45 trading days? clear;niter=1.0E5; % number of iterationsbelow=repmat(0,1,niter); % set up storagerandn( seed ,0);for i=1:niterr=normrnd(0.05/253,0.23/sqrt(253),1,45); % generate random numberslogPrice=log(1.0E6)+cumsum(r);minlogP=min(logPrice); %
5、 minmum price over next 45 daysbelow(i)=sum(minlogPlog(950000); endPro=mean(below)% P29 随机相位正弦波仿真% 1 time simulation w=2; N=1000; mu=2; sigma=3; s=rand(state); A=mu+sigma*randn(1,N); % A=normrnd(mu,sigma,1,N) theta=-pi+2*pi*rand(1,N); t=1:N; x=A.*cos(w*t+theta); capmu=mean(x) tao=1 x1=A.*cos(w*(t+ta
6、o)+theta); capgamma=mean(x-capmu).*(x1-capmu) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 17 页 - - - - - - - - - % m time simulation clear; w=2; N=1000; mu=2; sigma=3; m=500;capmu1=;capgamma1=; for i=1:m s=rand(state); A=mu+sigma*randn(1,N); theta=-pi+2*pi*r
7、and(1,N); t=1:N; x=A.*cos(w*t+theta); capmu=mean(x); capmu1=capmu1,capmu; tao=1; x1=A.*cos(w*(t+tao)+theta); capgamma=mean(x-capmu).*(x1-capmu); capgamma1=capgamma1,capgamma; end plot(1:m,capmu1,*,1:m,capgamma1,o) capmu=mean(capmu1); capgamma=mean(capgamma1); err1=mean(capmu1-0).2); gamma=(sigma2+mu
8、2)*cos(w*tao)/2; err2=mean(capgamma1-gamma).2); capmu,capgamma; err1, err2 % 输出： 0.0058 -2.7005 0.0065 0.0736 050100150200250300350400450500-4-3.5-3-2.5-2-1.5-1-0.500.5名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 17 页 - - - - - - - - - % P37 例 3.1.1 p1=poissc
9、df(5,10) p2=poisspdf(0,10) p1,p2 % 输出p1 =0.0671 p2 =4.5400e-005 ans =0.0671 0.0000 % P43 例 3.2.1 p3=poisspdf(9,12) % 输出p3 = 0.0874 % P43 例 3.2.2 p4=poisspdf(0,12) % 输出p4 = 6.1442e-006 % P39-40(Th3.1.1) Solve the difference equation system, find the solution % 输入：syms p0 p1 p2 ; S=dsolve(Dp0=-lamda*p
10、0,Dp1=-lamda*p1+lamda*p0,Dp2=-lamda*p2+lamda*p1, p0(0) = 1,p1(0) = 0,p2(0) = 0); S.p0,S.p1,S.p2 % 输出：ans = exp(-lamda*t), exp(-lamda*t)*t*lamda, 1/2*exp(-lamda*t)*t2*lamda2 % P43 泊松过程仿真% simulate 10 times clear; m=10; lamda=1; x=; for i=1:m s=exprnd(lamda,seed,1); x=x,exprnd(lamda); t1=cumsum(x); en
11、d x,t1 % 输出：名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 17 页 - - - - - - - - - ans = 0.6509 0.6509 2.4061 3.0570 0.1002 3.1572 0.1229 3.2800 0.8233 4.1033 0.2463 4.3496 1.9074 6.2570 0.4783 6.7353 1.3447 8.0800 0.8082 8.8882 % 输入：N=; for t=0:0.1:(t1(m)+1) if
12、 tt1(1) N=N,0; elseif tt1(2) N=N,1; elseif tt1(3) N=N,2; elseif tt1(4) N=N,3; elseif tt1(5) N=N,4; elseif tt1(6) N=N,5; elseif tt1(7) N=N,6; elseif tt1(8) N=N,7; elseif tt1(9) N=N,8; elseif tt1(10) N=N,9; else N=N,10; end end plot(0:0.1:(t1(m)+1),N,r-) % 输出：名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - -
13、 - - - - - - - 名师精心整理 - - - - - - - 第 5 页，共 17 页 - - - - - - - - - 012345678910012345678910% simulate 100 times clear; m=100; lamda=1; x=; for i=1:m s= rand(seed); x=x,exprnd(lamda); t1=cumsum(x); end x,t1 N=; for t=0:0.1:(t1(m)+1) if t=t1(i) & tt1(m) N=N,m; end end plot(0:0.1:(t1(m)+1),N,r-) % 输出：名
14、师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 17 页 - - - - - - - - - 01020304050607080901000102030405060708090100% P48 非齐次泊松过程仿真% simulate 10 timesclear; m=10; lamda=1; x=; for i=1:m s=rand( seed); % exprnd(lamda,seed,1); set seeds x=x,exprnd(lamda); t1=cumsum(
15、x); endx,t1 N=; T=; for t=0:0.1:(t1(m)+1) T=T,t.3; % time is adjusted, cumulative intensity function is t3. if t=t1(i) & tt1(m) N=N,m; endendplot(T,N,r-) % outputans = 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页，共 17 页 - - - - - - - - - 0.4220 0.4220 3.3323 3.
16、7543 0.1635 3.9178 0.0683 3.9861 0.3875 4.3736 0.2774 4.6510 0.2969 4.9479 0.9359 5.8838 0.4224 6.3062 1.7650 8.07120100200300400500600700800012345678910024681012x 105010203040506070809010010 times simulation 100 times simulation % P50 复合泊松过程仿真% simulate 100 times clear;niter=100; % iterate numberla
17、mda=1; % arriving ratet=input(Input a time:, s)for i=1:niter rand(state,sum(clock); x=exprnd(lamda); % interval time t1=x;while t1t x=x,exprnd(lamda); t1=sum(x); % arriving timeend t1=cumsum(x); y=trnd(4,1,length(t1); % rand(1,length(t1); gamrnd(1,1/2,1,length(t1); frnd(2,10,1,length(t1); t2=cumsum(
18、y);endx,t1,y,t2X=; m=length(t1);名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 8 页，共 17 页 - - - - - - - - - for t=0:0.1:(t1(m)+1)if t=t1(i) & tt1(m) X=X,t2(m);endendplot(0:0.1:(t1(m)+1),X,r-)0204060801001200510152025303540455002040608010012005101520253035404550跳跃度服从
19、0,1均匀分布情形跳跃度服从)2/1, 1(分布情形0102030405060708090-505101520跳跃度服从t （10）分布情形名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 9 页，共 17 页 - - - - - - - - - % Simulate the probability that sales revenue falls in some interval. (e.g. example 3.3.6 in teaching material) clear; ni
20、ter=1.0E4; % number of iterations lamda=6; % arriving rate (unit:minute) t=720; % 12 hours=720 minutes above=repmat(0,1,niter); % set up storage for i=1:niter rand(state,sum(clock); x=exprnd(lamda); % interval time n=1; while x=t n=n; else n=n+1; end end z=binornd(200,0.5,1,n); % generate n sales y=
21、sum(z); above(i)=sum(y432000); end pro=mean(above) Output: pro =0.3192 % Simulate the loss pro. For a Compound Poisson process clear; niter=1.0E3; % number of iterations lamda=1; % arriving rate t=input(Input a time:,s) below=repmat(0,1,niter); % set up storage for i=1:niter rand(state,sum(clock); x
22、=exprnd(lamda); % interval time n=1; while x=t 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 10 页，共 17 页 - - - - - - - - - n=n; else n=n+1; end end r=normrnd(0.05/253,0.23/sqrt(253),1,n); % generate n random jumps y=log(1.0E6)+cumsum(r); minX=min(y); % minmum return
23、 over next n jumps below(i)=sum(minXlog(950000); end pro=mean(below) Output: t=50, pro=0.45 % P86-87 (Example 5.1.5) markov chain chushivec0=0 0 1 0 0 0 P=0,1/2,1/2,0,0,0;1/2,0,1/2,0,0,0;1/4,1/4,0,1/4,1/4,0;0,0,1,0,0,0,;0,0,1/2,0,0,1/2;0,0,0,0,1,0 jueduivec1=chushivec0*P jueduivec2=chushivec0*(P2) %
24、 find 1 to n absolute-vector chushivec0=0 0 1 0 0 0; P=0,1/2,1/2,0,0,0;1/2,0,1/2,0,0,0;1/4,1/4,0,1/4,1/4,0;0,0,1,0,0,0,;0,0,1/2,0,0,1/2;0,0,0,0,1,0; n=10 t=1/6*ones(1 6); jueduivec=repmat(t,n 1); for k=1:n jueduiveck=chushivec0*(Pk); jueduivec(k,1:6)=jueduiveck end % Comparing two neighbour absolute
25、-vectors n=70; jueduivecn=chushivec0*(Pn) n=71; jueduivecn=chushivec0*(Pn) % Replace the first distribution, Comparing two neighbour absolute-vectors once more chushivec0=1/6 1/6 1/6 1/6 1/6 1/6; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 11 页，共 17 页 - - - - - -
26、- - - P=0,1/2,1/2,0,0,0;1/2,0,1/2,0,0,0;1/4,1/4,0,1/4,1/4,0;0,0,1,0,0,0,;0,0,1/2,0,0,1/2;0,0,0,0,1,0; n=70; jueduivecn=chushivec0*(Pn) n=71; jueduivecn=chushivec0*(Pn) % 赌博问题模拟（带吸收壁的随机游走：结束1次游走所花的时间及终止状态）a=5; p=1/2; m=0; while m0 & a0 & a15 m=m+1; r=2*binornd(1,p)-1; if r=-1 a=a-1; else a=a+1; enden
27、dif a=0 t1(1,n)=m; m1=m1+1; else t2(1,n)=m; m2=m2+1; endendfprintf(The average times of arriving 0 and 10 respectively are %d,%d.n,sum(t1,2)/m1,sum(t2,2)/m2); fprintf(The frequencies of arriving 0 and 10 respectively are %d,%d.n,m1/N, m2/N); % verify: fprintf(The probability of arriving 0 and its ap
28、proximate respectively are %d,%d.n, (p10*(1-p)5-p5*(1-p)10)/(p5*(p10-(1-p)10), m1/N); fprintf(The expectation of arriving 0 or 10 and its approximate respectively are %d,%d.n,5/(1-2*p)-10/(1-2*p)*(1-(1-p)5/p5)/(1-(1-p)10/p10), 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - -
29、 - - 第 13 页，共 17 页 - - - - - - - - - (sum(t1,2)+sum(t2,2)/N); 00.050.10.150.20.250.30.350.40.450.55101520253035404550pta0ta甲 的 预 期 输 光 时 间赌 博 平 均 持 续 时 间% 应用随机过程（04 年第一版） P125(example 5.29) 连续时间马尔可夫链PPT P29( 例 2) Solve the Kolmogorov difference equation,find the transation probabilities 输入：clear; sy
30、ms p00 p01 p10 p11 lamda mu; P=p00,p01;p10,p11; Q=-lamda,lamda;mu,-mu P*Q 输出：ans = -p00*lamda+p01*mu, p00*lamda-p01*mu -p10*lamda+p11*mu, p10*lamda-p11*mu 输入：p00,p01,p10,p11=dsolve(Dp00=-p00*lamda+p01*mu,Dp01=p00*lamda-p01*mu,Dp10=-p10*lamda+p11*mu,Dp11=p10*lamda-p11*mu,p00(0)=1,p01(0)=0,p10(0)=0,p1
31、1(0)=1) 输出：p00 = mu/(mu+lamda)+exp(-t*mu-t*lamda)*lamda/(mu+lamda) p01 = (lamda*mu/(mu+lamda)-exp(-t*mu-t*lamda)*lamda/(mu+lamda)*mu)/mu p10 = mu/(mu+lamda)-exp(-t*mu-t*lamda)*mu/(mu+lamda) p11 = (lamda*mu/(mu+lamda)+exp(-t*mu-t*lamda)*mu2/(mu+lamda)/mu 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - -
32、 - - - - - 名师精心整理 - - - - - - - 第 14 页，共 17 页 - - - - - - - - - % BPATH1 Brownian path simulation: for end randn( state,100) % set the state of randnT = 1; N = 500; dt = T/N; dW = zeros(1,N); % preallocate arrays .W = zeros(1,N); % for efficiencydW(1) = sqrt(dt)*randn; % first approximation outside
33、the loop .W(1) = dW(1); % since W(0) = 0 is not allowedfor j = 2:N dW(j) = sqrt(dt)*randn; % general increment W(j) = W(j-1) + dW(j); endplot(0:dt:T,0,W,r-) % plot W against txlabel(t, FontSize,16) ylabel(W(t), FontSize,16, Rotation,0) % BPATH2 Brownian path simulation: vectorized randn( state,100)
34、% set the state of randnT = 1; N = 500; dt = T/N; dW = sqrt(dt)*randn(1,N); % incrementsW = cumsum(dW); % cumulative sumplot(0:dt:T,0,W,r-) % plot W against txlabel(t, FontSize,16) ylabel(W(t), FontSize,16, Rotation,0) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1
35、5 页，共 17 页 - - - - - - - - - 00.10.20.30.40.50.60.70.80.91-0.500.51tW(t)%BPATH3 Function along a Brownian pathrandn( state,100) % set the state of randnT = 1; N = 500; dt = T/N; t = dt:dt:1; M = 1000; % M paths simultaneouslydW = sqrt(dt)*randn(M,N); % incrementsW = cumsum(dW,2); % cumulative sumU =
36、 exp(repmat(t,M 1) + 0.5*W); Umean = mean(U); plot(0,t,1,Umean,b-), hold on% plot mean over M pathsplot(0,t,ones(5,1),U(1:5,:),r-), hold off% plot 5 individual pathsxlabel(t, FontSize,16) ylabel(U(t), FontSize,16, Rotation,0, HorizontalAlignment, right) legend( mean of 1000 paths, 5 individual paths
37、,2) aveerr = norm(Umean - exp(9*t/8),inf) % sample error% 输出：名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 16 页，共 17 页 - - - - - - - - - 00.10.20.30.40.50.60.70.80.910.511.522.533.544.555.5tU(t)mean of 1000 paths5 individual pathsaveerr = 0.0504 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 17 页，共 17 页 - - - - - - - - -