gao_ins.m

％ PURPOSE : Demonstrate the differences between the following filters on the same problem: ％ ％ 1)

% PURPOSE : Demonstrate the differences between the following
% filters on an options pricing problem.
%           
%           1) Extended Kalman Filter  (EKF)
%           2) Unscented Kalman Filter (UKF)
%           3) Particle Filter         (PF)
%           4) PF with EKF proposal    (PFEKF)
%           5) 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.
P01 = 0.1;                  % EKF's initial variance of the
                            % interest rate.
P02 = 0.1;                  % EKF's initial variance of the volatility.

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_pfekf = 10*1e-5*eye(2);
R_pfekf = 1e-6*eye(2);

Q_pfukf = Q_pfekf;
R_pfukf = R_pfekf;
			    
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;
errorcEKF     = errorcTrivial;
errorpEKF     = errorcTrivial;
errorcUKF     = errorcTrivial;
errorpUKF     = errorcTrivial;
% LOAD THE DATA:
% =============
fprintf('\n')
fprintf('Loading the data')
fprintf('\n')
load c2925.prn;         load p2925.prn;
load c3025.prn;         load p3025.prn;
load c3125.prn;         load p3125.prn;
load c3225.prn;         load p3225.prn;
load c3325.prn;         load p3325.prn;
X=[2925; 3025; 3125; 3225; 3325];
[d1,i1]=sort(c2925(:,1));  Y1=c2925(i1,:);      Z1=p2925(i1,:);
[d2,i2]=sort(c3025(:,1));  Y2=c3025(i2,:);      Z2=p3025(i2,:);
[d3,i3]=sort(c3125(:,1));  Y3=c3125(i3,:);      Z3=p3125(i3,:);
[d4,i4]=sort(c3225(:,1));  Y4=c3225(i4,:);      Z4=p3225(i4,:);
[d5,i5]=sort(c3325(:,1));  Y5=c3325(i5,:);      Z5=p3325(i5,:);
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(3,:) = Y3(:,3)';   C(3,:) = Y3(:,2)';  P(3,:) = Z3(:,2)';
St(4,:) = Y4(:,3)';   C(4,:) = Y4(:,2)';  P(4,:) = Z4(:,2)';
St(5,:) = Y5(:,3)';   C(5,:) = Y5(:,2)';  P(5,:) = Z5(:,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:5
   Cox(i,:) = C(i,:) / X(i);
   Sox(i,:) = St(i,:) / X(i);
   Pox(i,:) = P(i,:) / X(i);
end
Cpred=zeros(T,5);
Ppred=zeros(T,5);

% 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 EKF and UKF ESTIMATION  %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%  ==============================  %%%%%%%%%%%%%%%%%%%%%

% INITIALISATION:
% ==============
mu_ekf = ones(2,T);       % EKF estimate of the mean of the states.
Inn = ones(2,2,T);        % Innovations Covariance.
Inn_ukf = Inn;
mu_ekf(1,1) = initr;
mu_ekf(2,1) = initsig;
P_ekf = ones(2,2,T);      % EKF estimate of the variance of the states.
for t=1:T
  P_ekf(:,:,t)= diag([P01 P02]);
end;
mu_ukf = mu_ekf;        % UKF estimate of the mean of the states.
P_ukf = P_ekf;          % UKF estimate of the variance of the states.

yPred = ones(2,T);      % One-step-ahead predicted values of y.
yPred_ukf = yPred;
mu_ekfPred = mu_ekf;    % EKF O-s-a estimate of the mean of the states.
PPred =eye(2);          % EKF O-s-a estimate of the variance of the states.
disp(' ');

for t=2:T,    
  fprintf('EKF & UKF : t = %i / %i  \r',t,T);
  fprintf('\n')
  
  % EKF PREDICTION STEP:
  % ==================== 
  mu_ekfPred(:,t) = feval('bsffun',mu_ekf(:,t-1),t);
  Jx = eye(2);  % Jacobian for bsffun.  
  PPred = Q + Jx*P_ekf(:,:,t-1)*Jx'; 
  
  % EKF CORRECTION STEP:
  % ====================
  yPred(:,t) = feval('bshfun',mu_ekfPred(:,t),u(:,t),t);

  % COMPUTE THE JACOBIAN:
  St  = u(1,t);            % Index price.
  tm  = u(2,t);            % Time to maturity.
  r   = mu_ekfPred(1,t);   % Risk free interest rate.
  sig = mu_ekfPred(2,t);   % Volatility.  
  d1 = (log(St) + (r+0.5*(sig^2))*tm ) / (sig * (tm^0.5));
  d2 = d1 - sig * (tm^0.5);  
  % Differentials of call price
  dcsig = St * sqrt(tm) * exp(-d1^2) / sqrt(2*pi);
  dcr   = tm * exp(-r*tm) * normcdf(d2);
  % Differentials of put price
  dpsig = dcsig;
  dpr   = -tm * exp(-r*tm) * normcdf(-d2);
  Jy = [dcr dpr; dcsig dpsig]'; % Jacobian for bshfun.

  % APPLY THE EKF UPDATE EQUATIONS:
  M = R + Jy*PPred*Jy';                 % Innovations covariance.
  Inn(:,:,t)=M;
  K = PPred*Jy'*inv(M);                 % Kalman gain.
  mu_ekf(:,t) = mu_ekfPred(:,t) + K*(y(:,t)-yPred(:,t));
  P_ekf(:,:,t) = PPred - K*Jy*PPred;
  
  % Full Unscented Kalman Filter step
  % =================================
  [mu_ukf(:,t),P_ukf(:,:,t),zab1,zab2,yPred_ukf(:,t),inov_ukf,Inn_ukf(:,:,t),K_ukf]=ukf(mu_ukf(:,t-1),P_ukf(:,:,t-1),u(:,t),Q_ukf,'ukf_bsffun',y(:,t),R_ukf,'ukf_bshfun',t,alpha,beta,kappa);
  
end;   % End of t loop.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-- CALCULATE PERFORMANCE --%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(3)
clf;
subplot(211)
p1=plot(1:T,mu_ekf(1,:),'r','linewidth',2);hold on;
p2=plot(1:T,mu_ukf(1,:),'b','linewidth',2);hold off;
legend([p1 p2],'ekf','ukf');
ylabel('Interest rate','fontsize',15)
subplot(212)
p1=plot(1:T,mu_ekf(2,:),'r','linewidth',2);hold on;
p2=plot(1:T,mu_ukf(2,:),'b','linewidth',2);hold off;
ylabel('Volatility','fontsize',15);
xlabel('Time (days)','fontsize',15)
zoom on;

% Transform innovations covariance for plotting.
Inn11=zeros(1,T);
Inn22=zeros(1,T);
Pekf11=zeros(1,T);
Pekf22=zeros(1,T);
for t=1:T,
  Inn11(t)=Inn(1,1,t);
  Inn22(t)=Inn(2,2,t);
  Inn11_ukf(t)=Inn_ukf(1,1,t);
  Inn22_ukf(t)=Inn_ukf(2,2,t);
  Pekf11(t)=P_ekf(1,1,t);
  Pekf22(t)=P_ekf(2,2,t);
  Pukf11(t)=P_ukf(1,1,t);
  Pukf22(t)=P_ukf(2,2,t);
end;

figure(1)
clf;
subplot(211)
plot(1:T,y(1,:),'r--',1:T,yPred(1,:),'b','linewidth',2);
hold on;
plot(1:T,yPred(1,:)+2*sqrt(Inn11),'k',1:T,yPred(1,:)-2*sqrt(Inn11),'k')
ylabel('Call price','fontsize',15)
legend('Actual price','Prediction');
axis([0 204 0.03 .22])
title('EKF');
subplot(212)
plot(1:T,y(2,:),'r--',1:T,yPred(2,:),'b','linewidth',2);
hold on;
plot(1:T,yPred(2,:)+2*sqrt(Inn22),'k',1:T,yPred(2,:)-2*sqrt(Inn22),'k')
ylabel('Put price','fontsize',15)
xlabel('Time (days)','fontsize',15)
axis([0 204 0 .06])
zoom on;
legend('Actual price','Prediction');

figure(2)
clf;
subplot(211)
plot(1:T,y(1,:),'r--',1:T,yPred_ukf(1,:),'b','linewidth',2);
hold on;
plot(1:T,yPred_ukf(1,:)+2*sqrt(Inn11_ukf),'k',1:T,yPred_ukf(1,:)-2*sqrt(Inn11_ukf),'k')
ylabel('Call price','fontsize',15)
legend('Actual price','Prediction');
axis([0 204 0.03 .22])
title('UKF');
subplot(212)
plot(1:T,y(2,:),'r--',1:T,yPred_ukf(2,:),'b','linewidth',2);
hold on;
plot(1:T,yPred_ukf(2,:)+2*sqrt(Inn22_ukf),'k',1:T,yPred_ukf(2,:)-2*sqrt(Inn22_ukf),'k')
ylabel('Put price','fontsize',15)
xlabel('Time (days)','fontsize',15)
axis([0 204 0 .06])
zoom on;
legend('Actual price','Prediction');
end