function NLL=NLLGP(Bet,Y,Bss);
%
%**************************************************************************
%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
%**************************************************************************
%
%**************************************************************************
%NEGATIVE LOG LIKELIHOOD FOR GENERALISED PARETO
%
%Bet     2p x 1 Spline parameter vector (for GP shape then GP scale)
%Y        n x 1 Responses
%Bss      n x p Spline basis evaluated at the covariate values corresponding to responses
%
%NLL      1 x 1 Negative log likelihood
%**************************************************************************

%Sizes
n=length(Y); %number of observations
nPrm=size(Bet,1)/2; %number of parameters

%Spline parameters and parameter values at responses
BetXi=Bet(1:nPrm); %Spline parameters for GP shape
BetSgm=Bet(nPrm+1:end); %Spline parameters for GP scale
Xi=Bss*BetXi; %GP shape at covariate values corresponding to responses
Sgm=Bss*BetSgm;  %GP scale at covariate values corresponding to responses
Nu=Sgm.*(1+Xi); %GP nu (in Nu = Sgm * (1 + Xi)), asymptotically uncorrelated with Xi 

%Tolerance
tol=1e-5; %tolerance value

%Initialise output
NLL_i=NaN(n,1); 

%Check validity conditions for likelihood
if any(Nu<eps & Xi>-1) || any(Nu>eps & Xi<-1) || any(Xi<-tol & (Y>-Nu./(Xi.*(1+Xi)))) || any(Xi<-0.5);
    NLL=1e10; %set likelihood to large "barrier" value
    return
end

%Calculate the negative log likelihood
I=abs(Xi)<tol;  %indicator for Gumbel case
nuI=Nu(I); %Nu for Gumbel case
nuNI=Nu(~I); %Nu for GP case
XiI=Xi(I); %Xi for Gumbel case
XiNI=Xi(~I); %Xi for GP case
NLL_i(I)=log(nuI)-log1p(XiI)+Y(I).*(1+XiI)./nuI; %Gumbel case
NLL_i(~I)=log(nuNI)-log1p(XiNI)+((1./XiNI) +1).*(log1p((XiNI./nuNI).*(1+XiNI).*Y(~I))); %GP case

%Calculate the sample penalised negative log likelihood
%Sum over individuals to get negative log likelihood
%Add roughness penalty terms for Xi and Sgm
NLL=sum(NLL_i);

return;