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 onylabel('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 benchmarkingfor 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.%Initializationfor 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]);  endend  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.endtime_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% meanserrorcTrivial_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);% varianceserrorcTrivial_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);