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(-1/(2
2、*(1-rho2)*(x-mu1)2/sigma12-2*rho*(x-mu1)*(y-mu2)/(sigma1*sigma2)+(y-mu2)2/sigma22);ezsurf(f)-6 -4 -2 0 2 4 600.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 deviation
3、of 0.23/year. These can be converted to rates per trading day by deviding by 253 and sqrt(253), respectively.Question 1: What is the probability that the value of the stock will be below $950,000 at the close day of at least one of the next 45 trading days? clear;niter=1.0E5; % number of iterationsb
4、elow=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); % minmum price over next 45 daysbelow(i)=sum(minlogP=t1(i) endendplot(0:0.1:(t1(m)+1),N,r-)% 输出：0 10 20 30 40 50 6
5、0 70 80 90 1000102030405060708090100% P48 非齐次泊松过程仿真% simulate 10 timesclear;m=10; lamda=1; x=; for i=1:ms=rand(seed); % exprnd(lamda,seed,1); set seedsx=x,exprnd(lamda);t1=cumsum(x);endx,t1N=; 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) endend
6、plot(T,N,r-)% outputans =0.4220 0.42203.3323 3.75430.1635 3.91780.0683 3.98610.3875 4.37360.2774 4.65100.2969 4.94790.9359 5.88380.4224 6.30621.7650 8.07120 100 200 300 400 500 600 700 8000123456789100 2 4 6 8 10 12x 105010203040506070809010010 times simulation 100 times simulation% P50 复合泊松过程仿真% si
7、mulate 100 timesclear;niter=100; % iterate numberlamda=1; % arriving ratet=input(Input a time:,s)for i=1:niterrand(state,sum(clock);x=exprnd(lamda); % interval timet1=x;while t1=t1(i) endendplot(0:0.1:(t1(m)+1),X,r-)0 20 40 60 80 100 120051015202530354045500 20 40 60 80 100 12005101520253035404550跳跃
8、度服从0,1均匀分布情形 跳跃度服从 分布情形)/,1(0 10 20 30 40 50 60 70 80 90-505101520跳跃度服从 t（10）分布情形% Simulate the probability that sales revenue falls in some interval. (e.g. example 3.3.6 in teaching material)clear; niter=1.0E4; % number of iterationslamda=6; % arriving rate (unit:minute)t=720; % 12 hours=720 minute
9、sabove=repmat(0,1,niter); % set up storagefor i=1:niterrand(state,sum(clock); x=exprnd(lamda); % interval timen=1;while x=t n=n;else n=n+1; end end z=binornd(200,0.5,1,n); % generate n salesy=sum(z);above(i)=sum(y432000); endpro=mean(above)Output: pro =0.3192% Simulate the loss pro. For a Compound P
10、oisson processclear; niter=1.0E3; % number of iterationslamda=1; % arriving ratet=input(Input a time:,s) below=repmat(0,1,niter); % set up storagefor i=1:niterrand(state,sum(clock);x=exprnd(lamda); % interval timen=1;while x=t n=n;else n=n+1; end endr=normrnd(0.05/253,0.23/sqrt(253),1,n); % generate
11、 n random jumpsy=log(1.0E6)+cumsum(r);minX=min(y); % minmum return over next n jumpsbelow(i)=sum(minX0 r=2*binornd(1,p)-1;if r=-1a=a-1;else a=a+1;endendif a=0t1(1,n)=m; m1=m1+1;elset2(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)
12、/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 approximate 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 respe
13、ctively are %d,%d.n,5/(1-2*p)-10/(1-2*p)*(1-(1-p)5/p5)/(1-(1-p)10/p10), (sum(t1,2)+sum(t2,2)/N);0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.55101520253035404550pta0 ta光光光光光光光光光光光光光光光光% 应用随机过程（04 年第一版）P125(example 5.29) 连续时间马尔可夫链 PPT P29(例 2)Solve the Kolmogorov difference equation,find the transati
14、on probabilities输入：clear;syms p00 p01 p10 p11 lamda mu; P=p00,p01;p10,p11;Q=-lamda,lamda;mu,-muP*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)
15、=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)/mup10 =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% BPATH1 Brownian path simulatio
16、n: forend 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 the loop .W(1) = dW(1); % since W(0) = 0 is not allowedfor j = 2:NdW(j) = sqrt(dt)*randn; % general
17、 incrementW(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) % set the state of randnT = 1; N = 500; dt = T/N;dW = sqrt(dt)*randn(1,N); % incrementsW = cumsum(dW); % c
18、umulative sumplot(0:dt:T,0,W,r-) % plot W against txlabel(t,FontSize,16)ylabel(W(t),FontSize,16,Rotation,0)0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-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 sim
19、ultaneouslydW = sqrt(dt)*randn(M,N); % incrementsW = cumsum(dW,2); % cumulative sumU = 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,2)aveerr = norm(Umean - exp(9*t/8),inf) % sample error% 输出：0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10.511.522.533.544.555.5tU(t)mean of 1000 paths5 individual pathsaveerr = 0.0504
