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  %%%%%%%%%%%%
%%%%%%  ======= EKF proposal ========  %%%%%%%% INITIALISATION:
xparticle_pfekf = 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_pfekf = cell(N,1);        % Particles for the covariance of x.
% Initialisation:
for i=1:N,
  xparticle_pfekf(1,1,i) = initr;   % sqrt(initr)*randn(1,1);
  xparticle_pfekf(2,1,i) = initsig; %sqrt(initsig)*randn(1,1);
  Pparticle_pfekf{i} = ones(2,2,T);
  for t=1:T,
    Pparticle_pfekf{i}(:,:,t)= diag([P01 P02]); 
  end;
end;  
xparticlePred_pfekf = ones(2,T,N);    % One-step-ahead predicted values of the states.
PparticlePred_pfekf = Pparticle_pfekf;    % One-step-ahead predicted values of P.
yPred_pfekf = ones(2,T,N);          % One-step-ahead predicted values of y.
w = ones(T,N);                      % Importance weights.
muPred_pfekf = ones(2,T);           % EKF O-s-a estimate of the mean of the states.
PPred_pfekf = ones(2,2);            % EKF O-s-a estimate of the variance of the states.
mu_pfekf = ones(2,T,N);             % EKF estimate of the mean of the states.
P_pfekf = ones(2,2,T);              % EKF estimate of the variance of the states.
disp(' ');
tic;                                % Initialize timer for benchmarking
for t=2:T,    
  fprintf('PF-EKF : t = %i / %i  \r',t,T);
  fprintf('\n')
  % PREDICTION STEP:
  % We use the EKF as proposal.
  for i=1:N,
    muPred_pfekf(:,t) = feval('bsffun',xparticle_pfekf(:,t-1,i),t);
    Jx = eye(2);                                 % Jacobian for ffun.
    PPred_pfekf = Q_pfekf + Jx*Pparticle_pfekf{i}(:,:,t-1)*Jx'; 
    yPredTmp = feval('bshfun',muPred_pfekf(:,t),u(:,t),t);
    % COMPUTE THE JACOBIAN:
    St  = u(1,t);              % Index price.
    tm  = u(2,t);              % Time to maturity.
    r   = muPred_pfekf(1,t);   % Risk free interest rate.
    sig = muPred_pfekf(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_pfekf + Jy*PPred_pfekf*Jy';                  % Innovations covariance.
   K = PPred_pfekf*Jy'*inv(M);                        % Kalman gain.
    mu_pfekf(:,t,i) = muPred_pfekf(:,t) + K*(y(:,t)-yPredTmp); % Mean of proposal.
    P_pfekf(:,:,t) = PPred_pfekf - K*Jy*PPred_pfekf;           % Variance of proposal.
    xparticlePred_pfekf(:,t,i) = mu_pfekf(:,t,i) + sqrtm(P_pfekf(:,:,t))*randn(2,1);
    PparticlePred_pfekf{i}(:,:,t) = P_pfekf(:,:,t);
  end;
  % EVALUATE IMPORTANCE WEIGHTS:
  % For our choice of proposal, the importance weights are give by:  
  for i=1:N,
    yPred_pfekf(:,t,i) = feval('bshfun',xparticlePred_pfekf(:,t,i),u(:,t),t);        
    lik = exp(-0.5*(y(:,t)-yPred_pfekf(:,t,i))'*inv(R)*(y(:,t)-yPred_pfekf(:,t,i)) ) + 1e-99;
    prior = exp(-0.5*(xparticlePred_pfekf(:,t,i)- xparticle_pfekf(:,t-1,i))'*inv(Q) * (xparticlePred_pfekf(:,t,i)-xparticle_pfekf(:,t-1,i) ))+ 1e-99;
    proposal = inv(sqrt(det(PparticlePred_pfekf{i}(:,:,t)))) * exp(-0.5*(xparticlePred_pfekf(:,t,i)-mu_pfekf(:,t,i))'*inv(PparticlePred_pfekf{i}(:,:,t)) * (xparticlePred_pfekf(:,t,i)-mu_pfekf(:,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_pfekf(:,t,:) = xparticlePred_pfekf(:,t,outIndex); % Keep particles with
                                                              % resampled indices.
  for i=1:N,
    Pparticle_pfekf{i} = PparticlePred_pfekf{outIndex(i)};  
  end;
end;   % End of t loop.
time_pfekf = toc;
%% 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:
yPFEKFmeanC=zeros(1,T);
yPFEKFmeanP=zeros(1,T);
for t=1:T,
  yPFEKFmeanC(t) = mean(yPred_pfekf(1,t,:));
  yPFEKFmeanP(t) = mean(yPred_pfekf(2,t,:));  
end;  
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));
errorcEKF(expr) =norm(C(104:204)-yPred(1,104:204));
errorpEKF(expr) =norm(P(104:204)-yPred(2,104:204));
errorcUKF(expr) =norm(C(104:204)-yPred_ukf(1,104:204));
errorpUKF(expr) =norm(P(104:204)-yPred_ukf(2,104:204));
errorcPF(expr) =norm(C(104:204)-yPFmeanC(104:204));
errorpPF(expr) =norm(P(104:204)-yPFmeanP(104:204));
errorcPFEKF(expr) =norm(C(104:204)-yPFEKFmeanC(104:204));
errorpPFEKF(expr) =norm(P(104:204)-yPFEKFmeanP(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(['EKF call       = ' num2str(errorcEKF(expr))]);
disp(['UKF call       = ' num2str(errorcUKF(expr))]);
disp(['PF call        = ' num2str(errorcPF(expr))]);
disp(['PF-EKF call    = ' num2str(errorcPFEKF(expr))]);
disp(['PF-UKF call    = ' num2str(errorcPFUKF(expr))]);
disp(['Trivial put    = ' num2str(errorpTrivial(expr))]);
disp(['EKF put        = ' num2str(errorpEKF(expr))]);
disp(['UKF put        = ' num2str(errorpUKF(expr))]);
disp(['PF put         = ' num2str(errorpPF(expr))]);
disp(['PF-EKF put     = ' num2str(errorpPFEKF(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);
p2=plot(bti:T,yPFEKFmeanC(bti:T),'r','linewidth',lw);
p3=plot(bti:T,yPFUKFmeanC(bti:T),'b','linewidth',lw); hold off;
ylabel('Call price','fontsize',15);
legend([p0 p1 p2 p3],'Actual price','PF prediction','PF-EKF 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);
p2=plot(bti:T,yPFEKFmeanP(bti:T),'r','linewidth',lw);
p3=plot(bti:T,yPFUKFmeanP(bti:T),'b','linewidth',lw); 
hold off;
ylabel('Put price','fontsize',15);
legend([p0 p1 p2 p3],'Actual price','PF prediction','PF-EKF 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);
errorcEKF_mean     = mean(errorcEKF);
errorcUKF_mean     = mean(errorcUKF);
errorcPF_mean      = mean(errorcPF);
errorcPFEKF_mean   = mean(errorcPFEKF);
errorcPFUKF_mean   = mean(errorcPFUKF);
errorpTrivial_mean = mean(errorpTrivial);
errorpEKF_mean     = mean(errorpEKF);
errorpUKF_mean     = mean(errorpUKF);
errorpPF_mean      = mean(errorpPF);
errorpPFEKF_mean   = mean(errorpPFEKF);
errorpPFUKF_mean   = mean(errorpPFUKF);
% variances
errorcTrivial_var = var(errorcTrivial);
errorcEKF_var     = var(errorcEKF);
errorcUKF_var     = var(errorcUKF);
errorcPF_var      = var(errorcPF);
errorcPFEKF_var   = var(errorcPFEKF);
errorcPFUKF_var   = var(errorcPFUKF);
errorpTrivial_var = var(errorpTrivial);
errorpEKF_var     = var(errorpEKF);
errorpUKF_var     = var(errorpUKF);
errorpPF_var      = var(errorpPF);
errorpPFEKF_var   = var(errorpPFEKF);
errorpPFUKF_var   = var(errorpPFUKF);
disp(' ');
disp('Mean and Variance of MSE ');
disp('-------------------------');
disp(' ');
disp(['Trivial call   : ' num2str(errorcTrivial_mean) ' (' num2str(errorcTrivial_var) ')']);
disp(['EKF call       : ' num2str(errorcEKF_mean) ' (' num2str(errorcEKF_var) ')']);
disp(['UKF call       : ' num2str(errorcUKF_mean) ' (' num2str(errorcUKF_var) ')']);
disp(['PF call        : ' num2str(errorcPF_mean) ' (' num2str(errorcPF_var) ')']);
disp(['PF-EKF call    : ' num2str(errorcPFEKF_mean) ' (' num2str(errorcPFEKF_var) ')']);
disp(['PF-UKF call    : ' num2str(errorcPFUKF_mean) ' (' num2str(errorcPFUKF_var) ')']);
disp(['Trivial put    : ' num2str(errorpTrivial_mean) ' (' num2str(errorpTrivial_var) ')']);
disp(['EKF put        : ' num2str(errorpEKF_mean) ' (' num2str(errorpEKF_var) ')']);
disp(['UKF put        : ' num2str(errorpUKF_mean) ' (' num2str(errorpUKF_var) ')']);
disp(['PF put         : ' num2str(errorpPF_mean) ' (' num2str(errorpPF_var) ')']);
disp(['PF-EKF put     : ' num2str(errorpPFEKF_mean) ' (' num2str(errorpPFEKF_var) ')']);
disp(['PF-UKF put     : ' num2str(errorpPFUKF_mean) ' (' num2str(errorpPFUKF_var) ')']);
figure(10);
subplot(211);
p1=semilogy(errorcTrivial,'k','linewidth',lw); hold on;
p2=semilogy(errorcEKF,'y','linewidth',lw);
p3=semilogy(errorcUKF,'g','linewidth',lw);
p4=semilogy(errorcPF,'m','linewidth',lw);
p5=semilogy(errorcPFEKF,'r','linewidth',lw);
p6=semilogy(errorcPFUKF,'b','linewidth',lw); hold off;
legend([p1 p2 p3 p4 p5 p6],'trivial','EKF','UKF','PF','PF-EKF','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;
p2=semilogy(errorpEKF,'y','linewidth',lw);
p3=semilogy(errorpUKF,'g','linewidth',lw);
p4=semilogy(errorpPF,'m','linewidth',lw);
p5=semilogy(errorpPFEKF,'r','linewidth',lw);
p6=semilogy(errorpPFUKF,'b','linewidth',lw); hold off;
legend([p1 p2 p3 p4 p5 p6],'trivial','EKF','UKF','PF','PF-EKF','PF-UKF');
ylabel('MSE','fontsize',12);
xlabel('experiment','fontsize',12);
title('PUT Options Mean Prediction Error','fontsize',14);