blackscholes.m

介绍了kf,ekf,ukf ,pf ,upf 的程序代码

% PURPOSE : Demonstrate the differences between the following
% filters on an options pricing problem.
%           
%           1) Particle Filter         (PF)
%           2)PF with UKF proposal    (PFUKF)


clear all;
echo off;
path('./ukf',path);
path('./data',path);


% INITIALISATION AND PARAMETERS:
% ==============================
doPlot = 0;                 % 1 plot online. 0 = only plot at the end.
g1 = 3;                     % Paramater of Gamma transition prior.
g2 = 2;                     % Parameter of Gamman transition prior.
                            % Thus mean = 3/2 and var = 3/4.
T = 204;                    % Number of time steps.
R = diag([1e-5 1e-5]);      % EKF's measurement noise variance. 
Q = diag([1e-7 1e-5]);      % EKF's process noise variance.

N = 10;                     % Number of particles.
optionNumber = 1;           % There are 5 pairs of options.
resamplingScheme = 1;       % The possible choices are
                            % systematic sampling (2),
                            % residual (1)
                            % and multinomial (3). 
                            % They're all O(N) algorithms. 

P01_ukf = 0.1
P02_ukf = 0.1

Q_ukf = Q;
R_ukf = R;

initr=.01
initsig=.15

Q_pfukf = 10*1e-5*eye(2);
R_pfukf = 1e-6*eye(2);
			    
alpha = 1;                     % UKF : point scaling parameter
beta  = 2;                     % UKF : scaling parameter for higher order terms of Taylor series expansion 
kappa = 1;                     % UKF : sigma point selection scaling parameter (best to leave this = 0)

no_of_experiments = 1;         % Number of times the experiment is
                               % repeated (for statistical purposes).

% DATA STRUCTURES FOR RESULTS
% ===========================

errorcTrivial = zeros(no_of_experiments,1);
errorpTrivial = errorcTrivial;
errorcPF      = errorcTrivial;
errorpPF      = errorcTrivial;
errorcPFUKF   = errorcTrivial;
errorpPFUKF   = errorcTrivial;


% LOAD THE DATA:
% =============
fprintf('\n')
fprintf('Loading the data')
fprintf('\n')
load c2925.prn;         load p2925.prn;
load c3025.prn;         load p3025.prn;
X=[2925; 3025];
[d1,i1]=sort(c2925(:,1));  Y1=c2925(i1,:);      Z1=p2925(i1,:);
[d2,i2]=sort(c3025(:,1));  Y2=c3025(i2,:);      Z2=p3025(i2,:);
d=Y1(:,1); 
% d - date to maturity.
St(1,:) = Y1(:,3)';   C(1,:) = Y1(:,2)';  P(1,:) = Z1(:,2)';
St(2,:) = Y2(:,3)';   C(2,:) = Y2(:,2)';  P(2,:) = Z2(:,2)';
% St - Stock price.
% C - Call option price.
% P - Put Option price.
% X - Strike price.
% Normalise with respect to the strike price:
for i=1:2
   Cox(i,:) = C(i,:) / X(i);
   Sox(i,:) = St(i,:) / X(i);
   Pox(i,:) = P(i,:) / X(i);
end
Cpred=zeros(T,2);
Ppred=zeros(T,2);

% PLOT THE LOADED DATA:
% ====================
figure(1)
clf;
plot(Cox');
ylabel('Call option prices','fontsize',15);
xlabel('Time to maturity','fontsize',15);
fprintf('\n')
fprintf('Press a key to continue')  
pause;
fprintf('\n')
fprintf('\n')
fprintf('Training has started')
fprintf('\n')

% SPECIFY THE INPUTS AND OUTPUTS
% ==============================
ii=optionNumber;   % Only one call price. Change 1 to 3, etc. for other prices.
X = X(ii,1);
St = Sox(ii,1:T);
C = Cox(ii,1:T);
P = Pox(ii,1:T);
counter=1:1:T;
tm = (224*ones(size(counter))-counter)/260;
u = [St' tm']';
y = [C' P']'; % Call and put prices.


% MAIN LOOP
% =========
for expr=1:no_of_experiments,

  rand('state',sum(100*clock));   % Shuffle the pack!
  randn('state',sum(100*clock));   % Shuffle the pack!  

%%%%%%%%%%%%%%%  PERFORM SEQUENTIAL MONTE CARLO  %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%  ==============================  %%%%%%%%%%%%%%%%%%%%%

% INITIALISATION:
% ==============
xparticle_pf = ones(2,T,N);      % These are the particles for the estimate
                                 % of x. Note that there's no need to store
                                 % them for all t. We're only doing this to
                                 % show you all the nice plots at the end.
xparticlePred_pf = ones(2,T,N);    % One-step-ahead predicted values of the states.
yPred_pf = ones(2,T,N);            % One-step-ahead predicted values of y.
w = ones(T,N);                   % Importance weights.
% Initialisation:
for i=1:N,
  xparticle_pf(1,1,i) = initr; % sqrt(initr)*randn(1,1);
  xparticle_pf(2,1,i) = initsig; %sqrt(initsig)*randn(1,1);
end;
disp(' ');
 
tic;                             % Initialize timer for benchmarking

for t=2:T,    
  fprintf('PF :  t = %i / %i  \r',t,T);
  fprintf('\n')
  
  % PREDICTION STEP:
  % ================ 
  % We use the transition prior as proposal.
  for i=1:N,
    xparticlePred_pf(:,t,i) = feval('bsffun',xparticle_pf(:,t-1,i),t) + sqrtm(Q)*randn(2,1);    
  end;

  % EVALUATE IMPORTANCE WEIGHTS:
  % ============================
  % For our choice of proposal, the importance weights are give by:  
  for i=1:N,
    yPred_pf(:,t,i) = feval('bshfun',xparticlePred_pf(:,t,i),u(:,t),t);        
    lik = exp(-0.5*(y(:,t)-yPred_pf(:,t,i))'*inv(R)*(y(:,t)-yPred_pf(:,t,i)) ) + 1e-99; % Deal with ill-conditioning.
    w(t,i) = lik;    
  end;  
  w(t,:) = w(t,:)./sum(w(t,:));                % Normalise the weights.
  
  % SELECTION STEP:
  % ===============
  % Here, we give you the choice to try three different types of
  % resampling algorithms. Note that the code for these algorithms
  % applies to any problem!
  if resamplingScheme == 1
    outIndex = residualR(1:N,w(t,:)');        % Residual resampling.
  elseif resamplingScheme == 2
    outIndex = systematicR(1:N,w(t,:)');      % Systematic resampling.
  else  
    outIndex = multinomialR(1:N,w(t,:)');     % Multinomial resampling.  
  end;
  xparticle_pf(:,t,:) = xparticlePred_pf(:,t,outIndex); % Keep particles with
                                                    % resampled indices.
end;   % End of t loop.

time_pf = toc;    % How long did this take?

% Compute posterior mean predictions:
yPFmeanC=zeros(1,T);
yPFmeanP=zeros(1,T);
for t=1:T,
  yPFmeanC(t) = mean(yPred_pf(1,t,:));
  yPFmeanP(t) = mean(yPred_pf(2,t,:));  
end;  





figure(4)
clf;
domain = zeros(T,1);
range = zeros(T,1);
thex=[0:1e-3:20e-3];
hold on
ylabel('Time (t)','fontsize',15)
xlabel('r_t','fontsize',15)
zlabel('p(r_t|S_t,t_m,C_t,P_t)','fontsize',15)
for t=11:20:200,
  [range,domain]=hist(xparticle_pf(1,t,:),thex);
  waterfall(domain,t,range/sum(range));
end;
view(-30,80);
rotate3d on;
a=get(gca);
set(gca,'ygrid','off');

figure(5)
clf;
domain = zeros(T,1);
range = zeros(T,1);
thex=[0.1:1e-2:0.25];
hold on
ylabel('Time (t)','fontsize',15)
xlabel('r_t','fontsize',15)
zlabel('p(\sigma_t|S_t,t_m,C_t,P_t)','fontsize',15)
%v=[0 1];
%caxis(v);
for t=11:20:200,
  [range,domain]=hist(xparticle_pf(2,t,:),thex);
  waterfall(domain,t,range/sum(range));
end;
view(-30,80);
rotate3d on;
a=get(gca);
set(gca,'ygrid','off');



%%%%%%%%%%%%%%%  PERFORM SEQUENTIAL MONTE CARLO  %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%  ======== UKF proposal ========  %%%%%%%%%%%%%%%%%%%%%

% INITIALISATION:
% ==============
xparticle_pfukf = ones(2,T,N);      % These are the particles for the estimate
                                    % of x. Note that there's no need to store
                                    % them for all t. We're only doing this to
                                    % show you all the nice plots at the end.
Pparticle_pfukf = cell(N,1);        % Particles for the covariance of x.

%Initialization
for i=1:N,
  xparticle_pfukf(1,1,i) = initr;   % sqrt(initr)*randn(1,1);
  xparticle_pfukf(2,1,i) = initsig; % sqrt(initsig)*randn(1,1);
  Pparticle_pfukf{i} = ones(2,2,T);
  for t=1:T,
    Pparticle_pfukf{i}(:,:,t) = diag([P01_ukf P02_ukf]);
  end
end  
xparticlePred_pfukf = ones(2,T,N);       % One-step-ahead predicted values of the states.
PparticlePred_pfukf = Pparticle_pfukf;   % One-step-ahead predicted values of P.
yPred_pfukf = ones(2,T,N);               % One-step-ahead predicted values of y.
w = ones(T,N);                           % Importance weights.
muPred_pfukf = ones(2,T);                % EKF O-s-a estimate of the mean of the states.
PPred_pfukf = ones(2,2);                 % EKF O-s-a estimate of the variance of the states.
mu_pfukf = ones(2,T,N);                  % EKF estimate of the mean of the states.
P_pfukf = ones(2,2,T);                   % EKF estimate of the variance of the states.

error=0;

disp(' ');

tic;

if (1),

for t=2:T,    
  fprintf('PF-UKF : t = %i / %i  \r',t,T);
  fprintf('\n')
  
  % PREDICTION STEP:
  % ================ 
  % We use the UKF as proposal.
  for i=1:N,
    % Call Unscented Kalman Filter
    [mu_pfukf(:,t,i),P_pfukf(:,:,t)]=ukf(xparticle_pfukf(:,t-1,i),Pparticle_pfukf{i}(:,:,t-1),u(:,t),Q_pfukf,'ukf_bsffun',y(:,t),R_pfukf,'ukf_bshfun',t,alpha,beta,kappa);
    xparticlePred_pfukf(:,t,i) = mu_pfukf(:,t,i) + sqrtm(P_pfukf(:,:,t))*randn(2,1);
    PparticlePred_pfukf{i}(:,:,t) = P_pfukf(:,:,t);
    
  end;

  % EVALUATE IMPORTANCE WEIGHTS:
  % ============================
  % For our choice of proposal, the importance weights are give by:  
  for i=1:N,
    yPred_pfukf(:,t,i) = feval('bshfun',xparticlePred_pfukf(:,t,i),u(:,t),t);

    lik = exp(-0.5*(y(:,t)-yPred_pfukf(:,t,i))'*inv(R)*(y(:,t)-yPred_pfukf(:,t,i)) ) + 1e-99;

    prior = exp(-0.5*(xparticlePred_pfukf(:,t,i)- xparticle_pfukf(:,t-1,i))'*inv(Q) * (xparticlePred_pfukf(:,t,i)-xparticle_pfukf(:,t-1,i) ))+ 1e-99;    
    proposal = inv(sqrt(det(PparticlePred_pfukf{i}(:,:,t)))) * exp(-0.5*(xparticlePred_pfukf(:,t,i)-mu_pfukf(:,t,i))'*inv(PparticlePred_pfukf{i}(:,:,t)) * (xparticlePred_pfukf(:,t,i)-mu_pfukf(:,t,i)))+ 1e-99;    
    
    w(t,i) = lik*prior/proposal;      
  end;  
  w(t,:) = w(t,:)./sum(w(t,:));                % Normalise the weights.
  
  % SELECTION STEP:
  % ===============
  % Here, we give you the choice to try three different types of
  % resampling algorithms. Note that the code for these algorithms
  % applies to any problem!
  if resamplingScheme == 1
    outIndex = residualR(1:N,w(t,:)');        % Residual resampling.
  elseif resamplingScheme == 2
    outIndex = systematicR(1:N,w(t,:)');      % Systematic resampling.
  else  
    outIndex = multinomialR(1:N,w(t,:)');     % Multinomial resampling.  
  end;
  xparticle_pfukf(:,t,:) = xparticlePred_pfukf(:,t,outIndex); % Keep particles with
                                              % resampled indices.
  for i=1:N,					      
    Pparticle_pfukf{i} = PparticlePred_pfukf{outIndex(i)};
  end
  
end;   % End of t loop.

end

time_pfukf = toc;


% Compute posterior mean predictions:
yPFUKFmeanC=zeros(1,T);
yPFUKFmeanP=zeros(1,T);
for t=1:T,
  yPFUKFmeanC(t) = mean(yPred_pfukf(1,t,:));
  yPFUKFmeanP(t) = mean(yPred_pfukf(2,t,:));  
end;  

errorcTrivial(expr) = norm(C(104:204)-C(103:203));
errorpTrivial(expr) = norm(P(104:204)-P(103:203));
errorcPF(expr) =norm(C(104:204)-yPFmeanC(104:204));
errorpPF(expr) =norm(P(104:204)-yPFmeanP(104:204));
errorcPFUKF(expr) =norm(C(104:204)-yPFUKFmeanC(104:204));
errorpPFUKF(expr) =norm(P(104:204)-yPFUKFmeanP(104:204));

disp(' ');
disp(['Experiment ' num2str(expr) ' of ' num2str(no_of_experiments) ' : Mean square errors sqrt(sum((errors).^2))']);
disp('------------------------------------------------------------');
disp(' ');
disp(['Trivial call   = ' num2str(errorcTrivial(expr))]);
disp(['PF call        = ' num2str(errorcPF(expr))]);
disp(['PF-UKF call    = ' num2str(errorcPFUKF(expr))]);
disp(['Trivial put    = ' num2str(errorpTrivial(expr))]);
disp(['PF put         = ' num2str(errorpPF(expr))]);
disp(['PF-UKF put     = ' num2str(errorpPFUKF(expr))]);


figure(9)
bti=20;
lw=2;
clf;
subplot(211)
p0=plot(bti:T,y(1,bti:T),'k-o','linewidth',lw); hold on;
p1=plot(bti:T,yPFmeanC(bti:T),'m','linewidth',lw);
p3=plot(bti:T,yPFUKFmeanC(bti:T),'b','linewidth',lw); hold off;
ylabel('Call price','fontsize',15);
legend([p0 p1 p3],'Actual price','PF prediction','PF-UKF prediction');
v=axis;
axis([bti T v(3) v(4)]);
subplot(212)
p0=plot(bti:T,y(2,bti:T),'k-o','linewidth',lw); hold on;
p1=plot(bti:T,yPFmeanP(bti:T),'m','linewidth',lw);
p3=plot(bti:T,yPFUKFmeanP(bti:T),'b','linewidth',lw); hold off;
ylabel('Put price','fontsize',15);
legend([p0 p1 p3],'Actual price','PF prediction','PF-UKF prediction');
xlabel('Time (days)','fontsize',15)
v=axis;
axis([bti T v(3) v(4)]);
zoom on;

end   % END OF MAIN LOOP

% CALCULATE MEAN AND VARIANCE OF EXPERIMENT RESULTS

% means
errorcTrivial_mean = mean(errorcTrivial);
errorcPF_mean      = mean(errorcPF);
errorcPFUKF_mean   = mean(errorcPFUKF);
errorpTrivial_mean = mean(errorpTrivial);
errorpPF_mean      = mean(errorpPF);
errorpPFUKF_mean   = mean(errorpPFUKF);

% variances
errorcTrivial_var = var(errorcTrivial);
errorcPF_var      = var(errorcPF);
errorcPFUKF_var   = var(errorcPFUKF);
errorpTrivial_var = var(errorpTrivial);
errorpPF_var      = var(errorpPF);
errorpPFUKF_var   = var(errorpPFUKF);

disp(' ');
disp('Mean and Variance of MSE ');
disp('-------------------------');
disp(' ');
disp(['Trivial call   : ' num2str(errorcTrivial_mean) ' (' num2str(errorcTrivial_var) ')']);
disp(['PF call        : ' num2str(errorcPF_mean) ' (' num2str(errorcPF_var) ')']);
disp(['PF-UKF call    : ' num2str(errorcPFUKF_mean) ' (' num2str(errorcPFUKF_var) ')']);
disp(['Trivial put    : ' num2str(errorpTrivial_mean) ' (' num2str(errorpTrivial_var) ')']);
disp(['PF put         : ' num2str(errorpPF_mean) ' (' num2str(errorpPF_var) ')']);
disp(['PF-UKF put     : ' num2str(errorpPFUKF_mean) ' (' num2str(errorpPFUKF_var) ')']);

figure(10);
subplot(211);
p1=semilogy(errorcTrivial,'k','linewidth',lw); hold on;
p4=semilogy(errorcPF,'m','linewidth',lw);
p6=semilogy(errorcPFUKF,'b','linewidth',lw); hold off;
legend([p1 p4 p6],'trivial','PF','PF-UKF');
ylabel('MSE','fontsize',12);
xlabel('experiment','fontsize',12);
title('CALL Options Mean Prediction Error','fontsize',14);
subplot(212);
p1=semilogy(errorpTrivial,'k','linewidth',lw); hold on;
p4=semilogy(errorpPF,'m','linewidth',lw);
p6=semilogy(errorpPFUKF,'b','linewidth',lw); hold off;
legend([p1 p4 p6],'trivial','PF','PF-UKF');
ylabel('MSE','fontsize',12);
xlabel('experiment','fontsize',12);
title('PUT Options Mean Prediction Error','fontsize',14);