随机过程matlab程序.docx
- 文档编号:23380597
- 上传时间:2023-05-16
- 格式:DOCX
- 页数:24
- 大小:304.65KB
随机过程matlab程序.docx
《随机过程matlab程序.docx》由会员分享,可在线阅读,更多相关《随机过程matlab程序.docx(24页珍藏版)》请在冰豆网上搜索。
随机过程matlab程序
%PPT例2一维正态密度与二维正态密度
symsxy;
s=1;t=2;
mu1=0;mu2=0;sigma1=sqrt((1+s^2));sigma2=sqrt((1+t^2));
x=-6:
0.1:
6;
f1=1/sqrt(2*pi*sigma1)*exp(-(x-mu1).^2/(2*sigma1^2));
f2=1/sqrt(2*pi*sigma2)*exp(-(x-mu2).^2/(2*sigma2^2));
plot(x,f1,'r-',x,f2,'k-.')
rho=(1+s*t)/(sigma1*sigma2);
f=1/(2*pi*sigma1*sigma2*sqrt(1-rho^2))*exp(-1/(2*(1-rho^2))*((x-mu1)^2/sigma1^2-2*rho*(x-mu1)*(y-mu2)/(sigma1*sigma2)+(y-mu2)^2/sigma2^2));
ezsurf(f)
%%Thedailylogreturnsonthestockhaveameanof0.05/yearandastandarddeviationof0.23/year.Thesecanbeconvertedtoratespertradingdaybydevidingby253andsqrt(253),respectively.
Question1:
Whatistheprobabilitythatthevalueofthestockwillbebelow$950,000atthe
closedayofatleastoneofthenext45tradingdays?
clear;
niter=1.0E5;%numberofiterations
below=repmat(0,1,niter);%setupstorage
randn('seed',0);
fori=1:
niter
r=normrnd(0.05/253,0.23/sqrt(253),1,45);%generaterandomnumbers
logPrice=log(1.0E6)+cumsum(r);
minlogP=min(logPrice);%minmumpriceovernext45days
below(i)=sum(minlogP end Pro=mean(below) %P29随机相位正弦波仿真 %1timesimulation 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+tao)+theta); capgamma=mean((x-capmu).*(x1-capmu)) %mtimesimulation clear; w=2;N=1000;mu=2;sigma=3; m=500;capmu1=[];capgamma1=[]; fori=1: m s=rand('state'); A=mu+sigma*randn(1,N); theta=-pi+2*pi*rand(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=(sigma^2+mu^2)*cos(w*tao)/2; err2=mean((capgamma1-gamma).^2); [capmu,capgamma;err1,err2] %输出: 0.0058-2.7005 0.00650.0736 %P37例3.1.1 p1=poisscdf(5,10) p2=poisspdf(0,10) [p1,p2] %输出 p1=0.0671 p2=4.5400e-005 ans=0.06710.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)Solvethedifferenceequationsystem,findthesolution %输入: symsp0p1p2; S=dsolve('Dp0=-lamda*p0','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)*t^2*lamda^2] %P43泊松过程仿真 %simulate10times clear; m=10;lamda=1;x=[]; fori=1: m s=exprnd(lamda,'seed',1); x=[x,exprnd(lamda)]; t1=cumsum(x); end [x',t1'] %输出: ans= 0.65090.6509 2.40613.0570 0.10023.1572 0.12293.2800 0.82334.1033 0.24634.3496 1.90746.2570 0.47836.7353 1.34478.0800 0.80828.8882 %输入: N=[]; fort=0: 0.1: (t1(m)+1) ift (1) N=[N,0]; elseift (2) N=[N,1]; elseift N=[N,2]; elseift N=[N,3]; elseift N=[N,4]; elseift N=[N,5]; elseift N=[N,6]; elseift N=[N,7]; elseift N=[N,8]; elseift N=[N,9]; else N=[N,10]; end end plot(0: 0.1: (t1(m)+1),N,'r-') %输出: %simulate100times clear; m=100;lamda=1;x=[]; fori=1: m s=rand('seed'); x=[x,exprnd(lamda)]; t1=cumsum(x); end [x',t1'] N=[]; fort=0: 0.1: (t1(m)+1) ift (1) N=[N,0]; end fori=1: (m-1) ift>=t1(i)&t N=[N,i]; end end ift>t1(m) N=[N,m]; end end plot(0: 0.1: (t1(m)+1),N,'r-') %输出: %P48非齐次泊松过程仿真 %simulate10times clear; m=10;lamda=1;x=[]; fori=1: m s=rand('seed');%exprnd(lamda,'seed',1);setseeds x=[x,exprnd(lamda)]; t1=cumsum(x); end [x',t1'] N=[];T=[]; fort=0: 0.1: (t1(m)+1) T=[T,t.^3];%timeisadjusted,cumulativeintensityfunctionist^3. ift (1) N=[N,0]; end fori=1: (m-1) ift>=t1(i)&t N=[N,i]; end end ift>t1(m) N=[N,m]; end end plot(T,N,'r-') %output ans= 0.42200.4220 3.33233.7543 0.16353.9178 0.06833.9861 0.38754.3736 0.27744.6510 0.29694.9479 0.93595.8838 0.42246.3062 1.76508.0712 10timessimulation100timessimulation %P50复合泊松过程仿真 %simulate100times clear; niter=100;%iteratenumber lamda=1;%arrivingrate t=input('Inputatime: ','s') fori=1: niter rand('state',sum(clock)); x=exprnd(lamda);%intervaltime t1=x; whilet1 x=[x,exprnd(lamda)]; t1=sum(x);%arrivingtime end 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(y); end [x',t1',y',t2'] X=[];m=length(t1); fort=0: 0.1: (t1(m)+1) ift (1) X=[X,0]; end fori=1: (m-1) ift>=t1(i)&t X=[X,t2(i)]; end end ift>t1(m) X=[X,t2(m)]; end end plot(0: 0.1: (t1(m)+1),X,'r-') 跳跃度服从[0,1]均匀分布情形跳跃度服从 分布情形 跳跃度服从t(10)分布情形 %%Simulatetheprobabilitythatsalesrevenuefallsinsomeinterval.(e.g.example3.3.6inteachingmaterial) clear; niter=1.0E4;%numberofiterations lamda=6;%arrivingrate(unit: minute) t=720;%12hours=720minutes above=repmat(0,1,niter);%setupstorage fori=1: niter rand('state',sum(clock)); x=exprnd(lamda);%intervaltime n=1; whilex x=x+exprnd(1/lamda);%arrivingtime ifx>=t n=n; else n=n+1; end end z=binornd(200,0.5,1,n);%generatensales y=sum(z); above(i)=sum(y>432000); end pro=mean(above) Output: pro=0.3192 %%Simulatethelosspro.ForaCompoundPoissonprocess clear; niter=1.0E3;%numberofiterations lamda=1;%arrivingrate t=input('Inputatime: ','s') below=repmat(0,1,niter);%setupstorage fori=1: niter rand('state',sum(clock)); x=exprnd(lamda);%intervaltime n=1; whilex x=x+exprnd(lamda);%arrivingtime ifx>=t n=n; else n=n+1; end end r=normrnd(0.05/253,0.23/sqrt(253),1,n);%generatenrandomjumps y=log(1.0E6)+cumsum(r); minX=min(y);%minmumreturnovernextnjumps below(i)=sum(minX end pro=mean(below) Output: t=50,pro=0.45 %P86-87(Example5.1.5)markovchain chushivec0=[001000] 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*(P^2) %find1tonabsolute-vector chushivec0=[001000]; 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([16]); jueduivec=repmat(t,[n1]); fork=1: n jueduiveck=chushivec0*(P^k); jueduivec(k,1: 6)=jueduiveck end %Comparingtwoneighbourabsolute-vectors n=70; jueduivecn=chushivec0*(P^n) n=71; jueduivecn=chushivec0*(P^n) %Replacethefirstdistribution,Comparingtwoneighbourabsolute-vectorsoncemore chushivec0=[1/61/61/61/61/61/6]; 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*(P^n) n=71; jueduivecn=chushivec0*(P^n) %赌博问题模拟(带吸收壁的随机游走: 结束1次游走所花的时间及终止状态) a=5;p=1/2; m=0; whilem<100 m=m+1; r=2*binornd(1,p)-1; ifr==-1 a=a-1; else a=a+1; end ifa==0|a==10 break; end end [ma] %赌博问题模拟(带吸收壁的随机游走: 结束N次游走所花的平均时间及终止状态分布规律) %p=q=1/2 p=1/2; m1=0;m2=0;N=1000; t1=0;t2=0; forn=1: 1: N m=0;a=5; whilea>0&a<10 m=m+1; r=2*binornd(1,p)-1; ifr==-1 a=a-1; else a=a+1; end end ifa==0 t1=t1+m;m1=m1+1; else t2=t2+m;m2=m2+1; end end fprintf('Theaveragetimesofarriving0and10respectivelyare%d,%d.\n',[t1/m1,t2/m2]); fprintf('Thefrequenciesofarriving0and10respectivelyare%d,%d.\n',[m1/N,m2/N]); %verify: fprintf('Theprobabilityofarriving0anditsapproximaterespectivelyare%d,%d.\n',[5/10,m1/N]); fprintf('Theexpectationofarriving0or10anditsapproximaterespectivelyare%d,%d.\n',[5*(10-5)/(2*p),(t1+t2)/N]); %p~=q p=1/4; m1=0;m2=0;N=1000; t1=zeros(1,N);t2=zeros(1,N); forn=1: 1: N m=0;a=5; whilea>0&a<15 m=m+1; r=2*binornd(1,p)-1; ifr==-1 a=a-1; else a=a+1; end end ifa==0 t1(1,n)=m;m1=m1+1; else t2(1,n)=m;m2=m2+1; end end fprintf('Theaveragetimesofarriving0and10respectivelyare%d,%d.\n',[sum(t1,2)/m1,sum(t2,2)/m2]); fprintf('Thefrequenciesofarriving0and10respectivelyare%d,%d.\n',[m1/N,m2/N]); %verify: fprintf('Theprobabilityofarriving0anditsapproximaterespectivelyare%d,%d.\n',[(p^10*(1-p)^5-p^5*(1-p)^10)/(p^5*(p^10-(1-p)^10)),m1/N]); fprintf('Theexpectationofarriving0or10anditsapproximaterespectivelyare%d,%d.\n',[5/(1-2*p)-10/(1-2*p)*(1-(1-p)^5/p^5)/(1-(1-p)^10/p^10),(sum(t1,2)+sum(t2,2))/N]); %应用随机过程(04年第一版)P125(example5.29)连续时间马尔可夫链PPTP29(例2) SolvetheKolmogorovdifferenceequation,findthetransationprobabilities 输入: clear; symsp00p01p10p11lamdamu; 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,p11(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)*mu^2/(mu+lamda))/mu %BPATH1Brownianpathsimulation: for…end randn('state',100)%setthestateofrandn T=1;N=500;dt=T/N; dW=zeros(1,N);%preallocatearrays... W=zeros(1,N);%forefficiency dW (1)=sqrt(dt)*randn;%firstapproximationoutsidetheloop... W (1)=dW (1);%sinceW(0)=0isnotallowed forj=2: N dW(j)=sqrt(dt)*randn;%gen
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 随机 过程 matlab 程序