随机过程matlab程序.doc

% PPT 例2 一维正态密度与二维正态密度
syms x y;
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)

% % The daily log returns on the stock have a mean of 0.05/year and a standard deviation 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 iterations
below=repmat(0,1,niter); % set up storage
randn('seed',0);
for i=1:niter
r=normrnd(0.05/253,0.23/sqrt(253),1,45); % generate random numbers
logPrice=log(1.0E6)+cumsum(r);
minlogP=min(logPrice); % minmum price over next 45 days
below(i)=sum(minlogP<log(950000));
end
Pro=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+tao)+theta);
capgamma=mean((x-capmu).*(x1-capmu))

% 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*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.0065 0.0736

% P37 例3.1.1
p1=poisscdf(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*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 泊松过程仿真
% 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);
end
[x',t1']

%输出：
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 t<t1(1)
N=[N,0];
elseif t<t1(2)
N=[N,1];
elseif t<t1(3)
N=[N,2];
elseif t<t1(4)
N=[N,3];
elseif t<t1(5)
N=[N,4];
elseif t<t1(6)
N=[N,5];
elseif t<t1(7)
N=[N,6];
elseif t<t1(8)
N=[N,7];
elseif t<t1(9)
N=[N,8];
elseif t<t1(10)
N=[N,9];
else
N=[N,10];
end
end
plot(0:0.1:(t1(m)+1),N,'r-')

%输出：

% 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(1)
N=[N,0];
end
for i=1:(m-1)
if t>=t1(i) & t<t1(i+1)
N=[N,i];
end
end
if t>t1(m)
N=[N,m];
end
end
plot(0:0.1:(t1(m)+1),N,'r-')

% 输出：

% P48 非齐次泊松过程仿真
% simulate 10 times
clear;
m=10; lamda=1; x=[];
for i=1:m
s=rand('seed'); % exprnd(lamda,'seed',1); set seeds
x=[x,exprnd(lamda)];
t1=cumsum(x);
end
[x',t1']

N=[]; T=[];
for t=0:0.1:(t1(m)+1)
T=[T,t.^3]; % time is adjusted, cumulative intensity function is t^3.
if t<t1(1)
N=[N,0];
end
for i=1:(m-1)
if t>=t1(i) & t<t1(i+1)
N=[N,i];
end
end
if t>t1(m)
N=[N,m];
end
end
plot(T,N,'r-')

% output
ans =
0.4220 0.4220
3.3323 3.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.0712

10 times simulation 100 times simulation

% P50 复合泊松过程仿真
% simulate 100 times
clear;
niter=100; % iterate number
lamda=1; % arriving rate
t=input('Input a time:','s')
for i=1:niter
rand('state',sum(clock));
x=exprnd(lamda); % interval time
t1=x;
while t1<t
x=[x,exprnd(lamda)];
t1=sum(x); % arriving time
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);
for t=0:0.1:(t1(m)+1)
if t<t1(1)
X=[X,0];
end
for i=1:(m-1)
if t>=t1(i) & t<t1(i+1)
X=[X,t2(i)];
end
end
if t>t1(m)
X=[X,t2(m)];
end
end
plot(0:0.1:(t1(m)+1),X,'r-')

跳跃度服从[0,1]均匀分布情形 跳跃度服从分布情形

跳跃度服从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 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
x=x+exprnd(1/lamda); % arriving time
if x>=t
n=n;
else
n=n+1;
end
end
z=binornd(200,0.5,1,n); % generate n sales
y=sum(z);
above(i)=sum(y>432000);
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=exprnd(lamda); % interval time
n=1;
while x<t
x=x+exprnd(lamda); % arriving time
if x>=t
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 over next n jumps
below(i)=sum(minX<log(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*(P^2)

% 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*(P^k);
jueduivec(k,1:6)=jueduiveck
end

% Comparing two neighbour absolute-vectors
n=70;
jueduivecn=chushivec0*(P^n)
n=71;
jueduivecn=chushivec0*(P^n)

% Replace the first distribution, Comparing two neighbour absolute-vectors once more
chushivec0=[1/6 1/6 1/6 1/6 1/6 1/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;
while m<100
m=m+1;
r=2*binornd(1,p)-1;
if r==-1
a=a-1;
else
a=a+1;
end
if a==0|a==10
break;
end
end
[m a]

% 赌博问题模拟（带吸收壁的随机游走：结束N次游走所花的平均时间及终止状态分布规律）
% p=q=1/2
p=1/2;
m1=0; m2=0; N=1000;
t1=0;t2=0;
for n=1:1:N
m=0; a=5;
while a>0 & a<10
m=m+1;
r=2*binornd(1,p)-1;
if r==-1
a=a-1;
else
a=a+1;
end
end
if a==0
t1=t1+m; m1=m1+1;
else
t2=t2+m; m2=m2+1;
end
end
fprintf('The average times of arriving 0 and 10 respectively are %d,%d.\n',[t1/m1,t2/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', [5/10, m1/N]);
fprintf('The expectation of arriving 0 or 10 and its approximate respectively are %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);
for n=1:1:N
m=0;a=5;
while a>0 & a<15
m=m+1;
r=2*binornd(1,p)-1;
if r==-1
a=a-1;
else
a=a+1;
end
end
if a==0
t1(1,n)=m; m1=m1+1;
else
t2(1,n)=m; m2=m2+1;
end
end
fprintf('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 approximate respectively are %d,%d.\n', [(p^10*(1-p)^5-p^5*(1-p)^10)/(p^5*(p^10-(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/p^5)/(1-(1-p)^10/p^10), (sum(t1,2)+sum(t2,2))/N]);

% 应用随机过程（04年第一版）P125(example 5.29) 连续时间马尔可夫链PPT P29(例2)
Solve the Kolmogorov difference equation,find the transation probabilities
输入：
clear;
syms 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,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

% BPATH1 Brownian path simulation: for…end
randn('state',100) % set the state of randn
T = 1; N = 500; dt = T/N;
dW = zeros(1,N); % preallocate arrays ...
W = zeros(1,N); % for efficiency

dW(1) = sqrt(dt)*randn; % first approximation outside the loop ...
W(1) = dW(1); % since W(0) = 0 is not allowed
for j = 2:N
dW(j) = sqrt(dt)*randn; % general increment
W(j) = W(j-1) + dW(j);
end

plot([0:dt:T],[0,W],'r-') % plot W against t

xlabel('t','FontSize',16)
ylabel('W(t)','FontSize',16,'Rotation',0)

% BPATH2 Brownian path simulation: vectorized

randn('state',100) % set the state of randn
T = 1; N = 500; dt = T/N;

dW = sqrt(dt)*randn(1,N); % increments
W = cumsum(dW); % cumulative sum

plot([0:dt:T],[0,W],'r-') % plot W against t

xlabel('t','FontSize',16)
ylabel('W(t)','FontSize',16,'Rotation',0)

%BPATH3 Function along a Brownian path

randn('state',100) % set the state of randn
T = 1; N = 500; dt = T/N; t = [dt:dt:1];

M = 1000; % M paths simultaneously
dW = sqrt(dt)*randn(M,N); % increments
W = cumsum(dW,2); % cumulative sum
U = exp(repmat(t,[M 1]) + 0.5*W);
Umean = mean(U);
plot([0,t],[1,Umean],'b-'), hold on % plot mean over M paths
plot([0,t],[ones(5,1),U(1:5,:)],'r--'), hold off % plot 5 individual paths

xlabel('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

% 输出：

aveerr = 0.0504