ekfdemo2.m

李子滤波器和神经网络在股票系统中的态势估计与系统盈亏趋势估计中的应用

% PURPOSE :  Sequential Tracking of Options Prices with an MLP model.
%       The model is estimated both with a plain extended Kalman filter 
%       (EKF) and an EKF where the process noise covariance is estimated 
%       adaptively.
             
% AUTHORS : Nando de Freitas and Mahesan Niranjan - Thanks for the acknowledgement :-)
% DATE    : 08-09-98

% LOAD THE DATA:
% =============

fprintf('\n')
fprintf('Loading the data')
fprintf('\n')
clear;
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);   N=length(d);
% d - date to maturity.
S(1,:) = Y1(:,3)';   C(1,:) = Y1(:,2)';  P(1,:) = Z1(:,2)';
S(2,:) = Y2(:,3)';   C(2,:) = Y2(:,2)';  P(2,:) = Z2(:,2)';
S(3,:) = Y3(:,3)';   C(3,:) = Y3(:,2)';  P(3,:) = Z3(:,2)';
S(4,:) = Y4(:,3)';   C(4,:) = Y4(:,2)';  P(4,:) = Z4(:,2)';
S(5,:) = Y5(:,3)';   C(5,:) = Y5(:,2)';  P(5,:) = Z5(:,2)';
% S - 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,:) = S(i,:) / X(i);
   Pox(i,:) = P(i,:) / X(i);
end
N = 204;
Cpred=zeros(N,5);
Ppred=zeros(N,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 the MLP with EKF')
fprintf('\n')

% SIMULATION:
% ==========

for ii=1:1   % Only one call price. Change 1 to 3, etc. for other prices.
  X = X(ii,1);
  S = Sox(ii,1:N);
  C = Cox(ii,1:N);
  P = Pox(ii,1:N);
  counter=1:1:N;
  tm = (224*ones(size(counter))-counter)/260;
  x = [S' tm']';

  % SIMULATION PARAMETERS:
  % =====================
  s1=6;             % Neurons in the hidden layer.
  s2=1;             % Neurons in the output layer.
  Q= 0.00001;       % Process noise covariance hyperparameter.
  R= 0.00001;       % Measurement noise covariance hyperparameter.
  initP = 10;       % Initial EKF covariance parameter.
  initVar = 1;      % Initial weights variance.
  window=10;        % Size of moving window to estimate Q.

  % PERFORM EKF ESTIMATION:
  % ======================
  tInit=clock;
  [p,theta,thetaR,PR,Innovations] = mlpekf(x,C,s1,s2,R,Q,initP,initVar,N);
  durationekf = etime(clock,tInit);

  % PERFORM EKF WITH EVIDENCE MAXIMISATION:
  % ======================================
  fprintf('\n')
  fprintf('Training the MLP with EKF/adaptive Q')
  fprintf('\n')
  tInit=clock;
  [p2,theta2,thetaR2,PRecord2,Innovations2,qplot] = mlpekfQ(x,C,s1,s2,R,Q,initP,initVar,window,N);
  durationekf2=etime(clock,tInit);

  errorT = norm(C(104:204)-p(104:204));
  errorT2 = norm(C(104:204)-p2(104:204));
  fprintf('\n')

  fprintf('Strike price = %d',X(ii))
  fprintf('\n')
  fprintf('EKF error    = %d',errorT)
  fprintf('\n')
  fprintf('EKF-Q erro r = %d',errorT2)
  fprintf('\n')
  fprintf('EKF duration   = %d seconds.\n',durationekf)
  fprintf('EKF-Q duration = %d seconds.\n',durationekf2)
  fprintf('\n')

  % PLOT FITTING:
  % ============

  figure(1)
  clf;
  subplot(221)
  plot(1:length(p),C,'g',1:length(p),p2,'b',1:length(p),p,'r')  
  legend('True value','EKF-Q estimate','EKF estimate');
  ylabel('One-step-ahead prediction','fontsize',15)
  xlabel('Time to maturity','fontsize',15)
  axis([0 204 0 .15]);
  subplot(222)
  plot(p,C,'r+',p2,C,'b+')
  ylabel('True value','fontsize',15)
  xlabel('Prediction','fontsize',15)
  legend('EKF estimate','EKF-Q estimate');
  hold on
  c=0.02:.01:.18;
  plot(c,c,'g');
  axis([0.02 .18 0.02 .18]);
  hold off
  subplot(223)
  plot(1:length(p),Innovations2)
  ylabel('Innovations variance for EKF-Q','fontsize',15)
  xlabel('Time','fontsize',15)
  subplot(224);
  plot(1:length(x),Innovations)
  ylabel('Innovations variance for EKF','fontsize',15)
  xlabel('Time','fontsize',15)

end %ii