function Rsl=KrnBS(Dat,Prm);
%**************************************************************************
%Illustration of non-stationary extreme value modelling in MATLAB accompanying  
%Threshold modeling of nonstationary extremes
%Paul Northrop, Philip Jonathan and David Randell
%in 
%Extreme Value Modeling and Risk Analysis: Methods and Applications
%Dipak Dey and Jun Yan (Editors)
%2015
%**************************************************************************
%
%**************************************************************************
%Kernel algorithm for bootstrapped MLE
%Dat          1 x 1 Dat structure
%Prm          1 x 1 Prm structure
%
%Rsl          1 x 1 Rsl structure
%**************************************************************************

OrgDat=Dat; %save a copy of the original data

%Bootstrapping
for iB=1:Prm.BSnRsm; %bootstrap loop
    
    if Prm.BSnRsm>1;
        fprintf('Bootstrap %g of %g\n',iB,Prm.BSnRsm);
    end
    
    %Resample the data
    if iB==1
        I=1:OrgDat.n; %first "bootstrap resample" is the original data
    else
        I=randi(OrgDat.n,OrgDat.n,1);
    end
    Dat.X=OrgDat.X(I);
    Dat.Y=OrgDat.Y(I);
    
    %Spline basis
    fprintf('\nMaking spline basis\n');
    [Bss,Pnl]=MakBss(Dat.X,Prm); %make spline basis
    nPrm=size(Pnl,1);  %number of model parameters;
    if Prm.BssIsPrd==0; %non-periodic
        Rsl.X=linspace(Prm.BssDmn(1),Prm.BssDmn(2),Prm.BssnDsp)'; %knot locations
    else; %periodic
        Rsl.X=linspace(Prm.BssDmn(1),Prm.BssDmn(2),Prm.BssnDsp+1)';
        Rsl.X=Rsl.X(2:end); %knot locations
    end;
    BssPrd=MakBss(Rsl.X,Prm); %basis for prediction
    
    %Quantile regression
    fprintf('\nEstimating quantile threshold\n');
    Bet0=ones(size(Pnl,1),1)*quantile(Dat.Y,Prm.QRNEP); %starting value just quantile of data;
    [Rsl.QR.OptBet,LOF]=CV(Prm.Rgh,@NLLLpl,Prm.CVnGrp,Dat.Y,Bss,Pnl,Bet0,1,Prm.QRNEP); %execute cross-validation loop
    if iB==1;
        Rsl.QR.LOF=LOF; %save predictive performance for actual sample
    end;
    %
    ThrDat=Bss*Rsl.QR.OptBet; %threshold at the data;
    if iB==1; 
        Rsl.QR.ThrDat=ThrDat;  %save threshold for actual sample
    end;
    Rsl.QR.MLE(:,iB)=BssPrd*Rsl.QR.OptBet; %MLE
 
    %Threshold exceedences
    IExc=Dat.Y>=ThrDat; %index of exceedences
    XExc=Dat.X(IExc);    

    %Poisson Rate
    fprintf('\nEstimating Poisson Rate of threshold exceedence\n')
    XCntEdg=linspace(Prm.BssDmn(1),Prm.BssDmn(2),Prm.PssnBin+1)';  %bin edges for possion counts;
    XCnt=(XCntEdg(1:end-1)+XCntEdg(2:end))/2;  %middle of bins for basis;
    Rsl.Pss.DltX=diff(XCntEdg(1:2));
    [BssCnt,~]=MakBss(XCnt,Prm); %make spline basis
    YCnt=histc(Dat.X(IExc),XCntEdg'); YCnt=YCnt(1:end-1); %count bin frequencies
    %
    Bet0=ones(size(Pnl,1),1)*mean(YCnt); %starting value mean count
    [Rsl.Pss.OptBet,LOF]=CV(Prm.Rgh,@NLLPss,Prm.CVnGrp,YCnt,BssCnt,Pnl,Bet0,1); %execute cross-validation loop
    if iB==1
        Rsl.Pss.LOF=LOF; %save predictive performance for actual sample
    end
    Rsl.Pss.MLE(:,iB)=BssPrd*Rsl.Pss.OptBet;; %MLE

    %Generalised Pareto
    fprintf('\nEstimating Generalised Pareto for threshold exceedence\n')
    YExc=Dat.Y(IExc)-ThrDat(IExc);
    [Bet0,LOFSSShp,LOFSSScl]=GPSS(Bss(IExc,:),Pnl,Prm.Rgh,YExc,XExc,Prm.GPHlfWdt,Prm.CVnGrp); %execute cross-validation loop
    %
    if iB==1; %save predictive performance for actual sample
        Rsl.GP.Shp.Str=BssPrd*Bet0(1:nPrm);
        Rsl.GP.Scl.Str=BssPrd*Bet0(nPrm+1:end);
        Rsl.GP.LOFSSShp=LOFSSShp;
        Rsl.GP.LOFSSScl=LOFSSScl;
    end;
    %
    [Rsl.GP.OptBet,LOF]=CV(Prm.Rgh,@NLLGP,Prm.CVnGrp,YExc,Bss(IExc,:),Pnl,Bet0,[Prm.GPSclFct,1]); %execute cross-validation loop
    Rsl.GP.Shp.MLE(:,iB)=BssPrd*Rsl.GP.OptBet(1:nPrm); %MLE
    Rsl.GP.Scl.MLE(:,iB)=BssPrd*Rsl.GP.OptBet(nPrm+1:end); %MLE
    %
    if iB==1;
        Rsl.GP.LOF=LOF; %save predictive performance for actual sample
    end;
    
end

%Figure with parameter estimates
figure(2);
clf;
subplot(2,2,1);
plot(OrgDat.X,OrgDat.Y,'k.','markersize',15);
hold on;
if size(Rsl.QR.MLE,2)==1; %only one realisation
    plot(Rsl.X,Rsl.QR.MLE,'r-','linewidth',2);
else; %multiple realisations, so plot quantiles
    plot(Rsl.X,quantile(Rsl.QR.MLE,0.025,2),'r--','linewidth',2);
    plot(Rsl.X,quantile(Rsl.QR.MLE,0.500,2),'r-','linewidth',2);
    plot(Rsl.X,quantile(Rsl.QR.MLE,0.975,2),'r--','linewidth',2);
end;
title('Threshold');
if isfield(OrgDat,'Thr')
    ThrTru=gpinv(Prm.QRNEP,OrgDat.Shp(OrgDat.X),OrgDat.Scl(OrgDat.X),OrgDat.Thr(OrgDat.X));
    plot(OrgDat.X,ThrTru,'g-','linewidth',2);
end
axis tight
%
subplot(2,2,2);
plot(XCnt,YCnt,'k.','markersize',15);
hold on;
if size(Rsl.Pss.MLE,2)==1; %only one realisation
    plot(Rsl.X,Rsl.Pss.MLE,'r-','linewidth',2);
else; %multiple realisations, so plot quantiles
    plot(Rsl.X,quantile(Rsl.Pss.MLE,0.025,2),'r--','linewidth',2);
    plot(Rsl.X,quantile(Rsl.Pss.MLE,0.500,2),'r-','linewidth',2);
    plot(Rsl.X,quantile(Rsl.Pss.MLE,0.975,2),'r--','linewidth',2);
end;
title('Rate');
%
hold on
if isfield(OrgDat,'Rat')
    plot(OrgDat.X,OrgDat.Rat(OrgDat.X)*Prm.QRNEP,'g-','linewidth',2);
end
axis tight
%
subplot(2,2,3);
cla; hold on;
if size(Rsl.QR.MLE,2)==1; %only one realisation
    plot(Rsl.X,Rsl.GP.Shp.MLE,'r-','linewidth',2);
else; %multiple realisations, so plot quantiles
    plot(Rsl.X,quantile(Rsl.GP.Shp.MLE,0.025,2),'r--','linewidth',2);
    plot(Rsl.X,quantile(Rsl.GP.Shp.MLE,0.500,2),'r-','linewidth',2);
    plot(Rsl.X,quantile(Rsl.GP.Shp.MLE,0.975,2),'r--','linewidth',2);
end;
title('Shape');
plot(Rsl.X,Rsl.GP.Shp.Str,'b--','linewidth',2);
if isfield(OrgDat,'Shp')
    plot(OrgDat.X,OrgDat.Shp(OrgDat.X),'g-','linewidth',2);
end
axis tight
%
subplot(2,2,4);
cla; hold on;
if size(Rsl.QR.MLE,2)==1; %only one realisation
    plot(Rsl.X,Rsl.GP.Scl.MLE,'r-','linewidth',2);
else; %multiple realisations, so plot quantiles
    plot(Rsl.X,quantile(Rsl.GP.Scl.MLE,0.025,2),'r--','linewidth',2);
    plot(Rsl.X,quantile(Rsl.GP.Scl.MLE,0.500,2),'r-','linewidth',2);
    plot(Rsl.X,quantile(Rsl.GP.Scl.MLE,0.975,2),'r--','linewidth',2);
end;
title('Scale');
plot(Rsl.X,Rsl.GP.Scl.Str,'b--','linewidth',2);
if isfield(OrgDat,'Scl')
    tScl=OrgDat.Scl(OrgDat.X)+OrgDat.Shp(OrgDat.X).*(Rsl.QR.ThrDat-OrgDat.Thr(OrgDat.X));
    plot(OrgDat.X,tScl,'g-','linewidth',2);
end
axis tight

%Figure with lack of fit
figure(3);
clf
subplot(2,2,1)
Rsl.QR.LOF(Rsl.QR.LOF>1e9)=NaN;
plot(Prm.Rgh,Rsl.QR.LOF,'k.-','linewidth',2)
set(gca,'xscale','log');
ylabel('QR')
%
subplot(2,2,2)
Rsl.Pss.LOF(Rsl.Pss.LOF>1e9)=NaN;
plot(Prm.Rgh,Rsl.Pss.LOF,'k.-','linewidth',2)
set(gca,'xscale','log');
ylabel('Pss')
%
subplot(4,2,5)
Rsl.GP.LOFSSShp(Rsl.GP.LOFSSShp>1e9)=NaN;
plot(Prm.Rgh, Rsl.GP.LOFSSShp,'k.-','linewidth',2)
set(gca,'xscale','log');
ylabel('StrShp')
%
subplot(4,2,7)
Rsl.GP.LOFSSScl(Rsl.GP.LOFSSScl>1e9)=NaN;
plot(Prm.Rgh,Rsl.GP.LOFSSScl,'k.-','linewidth',2)
set(gca,'xscale','log');
ylabel({'LOF';'StrScl'})
xlabel('Roughness penalty');
%
subplot(2,2,4)
Rsl.GP.LOF(Rsl.GP.LOF>1e9)=NaN;
plot(Prm.Rgh,Rsl.GP.LOF,'k.-','linewidth',2)
set(gca,'xscale','log');
ylabel('GP')

%Return value simulation
Rsl=RV(Rsl,Prm,Dat);

return;