% KernelFilterMAE
% y(t) = x(t) + epsilon(t)
% x(t) = x(t-1) + eta(t)
% N = # of particles. 
clear
bugtrap = 1
Reps = 1000  % Reps <= 1000
T = 100
N = 100   % N = 100, 1000, 10000, 100000
alpha = 1.1    % = 1.1, 1.3, 1.5, 1.7, 1.9
ceps = 1
ceta = .25
betaeps = .3
betaeta = 0
param = 1 
%  = 0 for continuous zeta parameterization
%  = 1 for mean-focus mu parameterization
seed = 1
rng(seed)
% Generate and plot random data
% Uses Nolan RNG etc.  Defaults to S0 parameterization (continuous) ,
%   but using Nolan S1(mean-focus) parm. 
tic
eps = stablernd(1000*T, [alpha betaeps ceps 0], param);
eps = reshape(eps, 1000, T);
eps = eps(1:Reps, :);
'eps done'
eta = stablernd(1000*T, [alpha betaeta ceta 0], param);
eta = reshape(eta, 1000, T);
eta = eta(1:Reps, :);
'eta done'
x = cumsum(eta')';
y = x + eps;    % x, y ~ Reps x T  Reps must be <= 1000
'x y done'
toc
% x, y ~ Reps x T

tic
ESSbar = nan(1,Reps);
ESSmin = nan(1,Reps);
MAE = nan(1,Reps);
plots = -1

% plots = 1
for r = 1:Reps 
  r = r      
  % use r as seed
   [muF, sigF, pF, logl, pci, ci, ESS, tictoc] =  ...
      StableKernelFilter(y(r,:) , alpha, ceps, ceta, betaeps, N, param, plots, r);
%   ESS ~ 1xT
%   muF etc ~ (T, N)
  ESSbar(r) = mean(ESS);
  ESSmin(r) = min(ESS);
  EF = dot(muF, pF, 2); % ~ Tx1  equivalently, EF = dot(muF', pF')'
  E = EF' - x(r,:); 
  MAE(r) = mean(abs(E(11:T)));
  % end % if r == 273 etc
%   if Reps < 11
%     MAEr = MAE(r)
%     figure
%     plot(1:T, [E; abs(E); MAEr*ones(1,T); zeros(1,T)])
%     title (['Errors, r = ' num2str(r)])
%     figure
%     plot(1:T, ESS)
%     title (['ESS, r = ' num2str(r)])
%   end % if Reps    
end % for r
'Reps done'
toc

tic
figure
plot(1:Reps, ESSbar, 'linewidth', 2)
title('ESSbar')
xlabel('replication')
grid on
figure
plot(1:Reps, ESSmin, 'linewidth', 2)
title('ESSmin')
xlabel('replication')
grid on
figure
plot(1:Reps, MAE, 'linewidth', 2)
title('MAE')
xlabel('replication')
grid on
'plots done'
toc
alpha = alpha
N=N
Reps=Reps


avgMAE = mean(MAE)
medMAE = median(MAE)
Q90MAE = quantile(MAE, .90)
maxMAE = max(MAE)
avgESSmin = mean(ESSmin)
minESSmin = min(ESSmin)
Q10ESSmin = quantile(ESSmin,.10)
medESSmin = median(ESSmin)

return


