% New Particle Filter Demo 1520351
% y(t) = x(t) + epsilon(t)
% x(t) = x(t-1) + eta(t)
% N = # of particles.  Must be even for Mixed Particulate Adaptive
% equilibration
clear

N = 100 
T = 100
alpha = 1.7
ceps = 1
ceta = .25
betaeps = .3
betaeta = 0
param = 1 
%  = 0 for continuous zeta parameterization
%  = 1 for mean-focus mu parameterization
plots = 1
% plots = 2 shows particles etc.
% plots = 1 shows data, filter, CIs, no particles

seed = 1520351 
rng(seed)
% Generate and plot random data
% Uses Nolan RNG etc.  Defaults to 0 parameterization (continuous)
eps = stablernd(T, [alpha betaeps ceps 0], param); 
eta = stablernd(T, [alpha betaeta ceta 0], param);
x = cumsum(eta); % x(1) = eta(1)
y = x + eps;
if plots > 0
  figure 
  plot(1:T, eta, 'bd', 1:T, eps, 'r*', 1:T, zeros(1, T), 'k-', 'markersize', 10)
  legend('eta', 'epsilon')
  grid on
  figure
  plot(1:T, x, 'b-d', 1:T, y, 'r*', 'linewidth', 1, 'markersize', 10)
  grid on
  legend ('x', 'y')
  xlabel('Time t')
end % if plots

seed = 1
rng(seed)
[xF, pF, logl, pci, ci, ESS, tictoc] = StableParticleFilter(y, alpha, ... 
  ceps, ceta, betaeps, N, param, plots, seed);
tictoc = tictoc





