function linmarkjumpdiff06fig1 % Fig. 5.1 Book Illustration for Linear Distributed-Jump Diffusion (3/2007) % SDE RNG Simulation with variable coefficients for t in [0,1] % with sample variation: % DX(t) = X(t)*(mu(t)*Dt + sig(t)*DW(t) + nu(Q)*DP(t), % X(0) = x0. % Or log-state: % DY(t) = (mu(t)-sig^2(t)/2)*Dt + sig(t)*DW(t) + Q*DP(t), % Y(0) = log(x0) and Q = ln(1+nu(Q)). % Generation is by summing Wiener increments DW of even spacing Dt % with Poisson jump increment added at correct time increment. % Sufficiently SMALL increments assumed, so zero-one jump law is % appropriate. % For demonstration purposes, Q will be assumed to be % (qdist =1) UNIFORMLY distributed on (qparm1,qparm2)=(a,b) % OR % (qdist=2) NORMALLY distributed with (qparm1,qparm2)=(muj,sj2). % Allows Separate Driver Input and Special Jump % or Diffusion Handling. clc % clear variables, but must come before globals, % else clears globals too. clf % clear figures fprintf('\nfunction linjumpdiff06fig1 OutPut:'); %%% Initialize input to jdsimulator with sample parameters: N = 1000; t0 = 0; T = 2.0; % Set initial time grid: Fixed Delta{t}. idiff = 1; ijump = 1; x0 = 1.0; qdist = 1; a = -2; b = +1; qparm1 = a; qparm2 = b; %e.g., Uniform %OR E.G., Normal distribution: %qdist = 2; muj = 0.28; sj2 = +0.15; qparm1 = muj; qparm2 = sj2; % set constant parameters. fprintf('\n N=%i; x0=%6.3f; t0=%6.3f; T=%6.3f;',N,x0,t0,T); fprintf('\n qdist=%i*; qparm1=%6.3f; qparm2=%6.3f;'... ,qdist,qparm1,qparm2); fprintf('\n * qdist=1 for uniform Q-distribution.'); fprintf('\n * qdist=2 for normal Q-distribution.'); % jdsimulator(idiff,ijump,qdist,qparm1,qparm2,N,x0,t0,T); % % END INPUT FOR JUMP-DIFFUSION SIMULATOR. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function jdsimulator(idiff,ijump,qdist,qparm1,qparm2,N,x0,t0,T) kfig = 0; % figure counter. dt = (T-t0)/N; % Get number of intervals/samples and time step. kjd = 4 - 2*idiff - ijump; NP = N + 1; % Number of plot points = number of time steps + 1. sqrtdt = sqrt(dt); % Set standard Wiener increment moments. tv = t0:dt:T; % Compute time vector; sv = zeros(size(tv)); ldtv = zeros(size(tv)); muv = mu(tv); % Get time-dependent coefficient vectors if idiff == 1, sv = sigma(tv); end if ijump == 1, ldtv = dt*lambda(tv); end muddt = (muv - sv.^2/2)*dt; % Get diffusion corrected drift term. if qdist == 1 % Average nu(Q)=exp(Q)-1 for UNIFORM Q-Dist. numean = (exp(qparm2)-exp(qparm1))/(qparm2-qparm1)-1; elseif qdist == 2 % Average nu(Q)=exp(Q)-1 for NORMAL Q-Dist. numean = exp(qparm1-qparm2/2)-1; end % Compute Theoretical Mean State Path % E[X(t+dt)] = X(t)*exp(E[dX(t)|X(t)=x]/x), x0 > 0: XM = zeros(1,NP); % preallocate mean state. XM(1) = x0; for i = 1:N XM(i+1) = XM(i)*exp(muv(i)*dt+numean*ldtv(i)); end kstates = 4; kv = [1,5,9,10]; % selected random states. XS = zeros(NP,kstates); % preallocate global state array. % Begin Sample Path Calculation: for k = 1:kstates % Test Multiple Simulated Sample Paths: if idiff == 1 randn('state',kv(k)); % Set initial normal state % for repeatability. DW = sqrtdt*randn(1,N); % Generate normal random vector % of N samples for DW(t). end if ijump == 1 rand('state',kv(k)); % Set initial uniform state % for repeatability. DU = rand(1,N); % Generate Uniform random vector DP(t) if qdist == 1 %Generate Uniform random mark vector Q samples. Q = qparm1+(qparm2-qparm1)*rand(1,N); elseif qdist == 2 %Generate Normal random mark vector Q samples. sj = sqrt(qparm2); Q = qparm1+sj*randn(1,N); end ul = (1-ldtv)/2; ur = 1-ul; % Set vector centered jump % probability thresholds, end YS = zeros(1,N+1); % preallocate state exponent for efficiency. XS(1,k) = x0; % Set kth initial state. for i = 1:N % Simulated Sample paths by Increment Accumulation: YS(i+1) = YS(i) + muddt(i); % Add dY-drift: % Add diffusion increment: if idiff == 1, YS(i+1) = YS(i+1)+ sv(i)*DW(i); end % Using centered part of uniform distribution, with % acceptance-rejection, to avoid open end point bias: if ijump == 1 if DU(i) <= ur(i) && DU(i) >= ul(i) % Jump if in [ul,ur] YS(i+1) = YS(i+1) + Q(i); % If jump, +Y-jump amplitude. end % Else no jump, so do not add anything. end XS(i+1,k) = x0*exp(YS(i+1));% Invert exponent to get state. end % i end % k % Sample Mean State: XSM = zeros(1,NP); for i = 1:NP XSM(i) = mean(XS(i,:)); end % Begin Plot: scrsize = get(0,'ScreenSize'); ss = 5.2; dss = 0.2; ssmin = 3.0; kfig = kfig + 1; stitle = {'Linear Mark-Jump-Diffusion Simulations' ... ,'Linear Diffusion Simulations' ... ,'Linear Mark-Jump Simulations'}; sylabel = {'X(t), Jump-Diffusion State','X(t), Diffusion State' ... ,'X(t), Jump State'}; slegend = {'X(t), State 1','X(t), State 5' ... ,'X(t), State 9','X(t), State 10'... ,'XM(t), th. Mean=E[X(t)]','XSM(t), Sample Mean'}; fprintf('\n\nFigure(%i): Linear Jump-Diffusion Simulations\n',kfig) figure(kfig) plot(tv,XS(1:NP,1),'k+-' ... ,tv,XS(1:NP,2),'k:' ... ,tv,XS(1:NP,3),'k--' ... ,tv,XS(1:NP,4),'k-.' ... ,tv,XM(1:NP),'k-',tv,XSM(1:NP),'b.-' ... ,'LineWidth',2); % Add for more States? title(stitle(kjd),'FontWeight','Bold','Fontsize',44); ylabel(sylabel(kjd),'FontWeight','Bold','Fontsize',44); xlabel('t, Time','FontWeight','Bold','Fontsize',44); hlegend=legend(slegend,'Location','NorthEast'); set(hlegend,'Fontsize',36,'FontWeight','Bold'); set(gca,'Fontsize',36,'FontWeight','Bold','linewidth',3); ss = max(ss - dss,ssmin); set(gcf,'Color','White','Position' ... ,[scrsize(3)/ss 60 scrsize(3)*0.60 scrsize(4)*0.80]); %[l,b,w,h] % % End JDSimulator Code % % linear Time-Dependent SDE Coefficient Functions: % (Change with application; functions must be vectorizable, % using vector element dot operations or vector functions.) % function v = mu(t) % drift coefficient example, change with applications: v = 0.1*sin(t); % end mu(t) % function v = sigma(t) % drift coefficient example, change with applications: v = 1.5*exp(-0.01*t); % end sigma(t) % function v = lambda(t) % drift coefficient example, change with applications: v = 3.0*exp(-t.*t); % end lambda(t)