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)